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ABSTRACT 

There has been a continuing problem of estimating a 
radar signal when the noise and the signal have the same 
power spectra. This is particularly troublesome when one 
tries to resolve two close targets with a tracking radar. 

The purpose of this research is to show theoretically that 
there are practical ways to solve the problem. 

Two approaches are introduced here. The first is to 
add a process filter to the radar, and the second is to 
design a new measuring technique with processing such that 
the signals will be separable. These two main approaches 
have led to several new or extended theories which were 
developed in the course of this work: 

The Probability filter. 

A modified Kalman filter (MKF) . 

A measurement modification technique for monopulse 
radar. 

A new approach to the design of a monopulse tracking 
radar . 

Simulations were performed to check the two major 
theories (the M.K.F. and the probability filter). According 
to the simulations we can say that a tracking radar can be 
modified in order to solve the problem of separating 
unresolved targets. 
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I. 



INTRODUCTION 



A. SEPARATION OF SIGNALS 

In many cases we would like to separate two signals 
from a measurement which is a function of their combination, 
such as their sura. One of the signals is usually "noise" 
and the other the "signal" or "information." Techniques of 
separating two signals if they have different power spectra 
are well known, for example, the Kalman filter or the Weiner 
filter. The case when the two signals have almost the same 
power spectra does not appear in the literature, mostly 
because each problem must be treated in a different way. 

The following are two basic approaches to this problem: 

(i) Use all of the statistical information about 
the signal to effect the separation. 

(ii) Design measuring techniques such that the two 
signals will be separable in some sense. 

The specific application to be addressed in this research is 
the tracking radar case where it is necessary to separate a 
target return from chaff. The sections which follow present 
a summary of the basic theory with reference to the litera- 
ture. Some new results have been achieved and are presented 
in the body of the text. They include; 

a) The separation of the two signals by means of a 
probability filter - Appendix D 

b) A modified Kalmak filter approach - Appendix A. 
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c) A new measurement modification technique for mono- 
pulse radar - Section (IV. A) . 

d) A new approach to a tracking radar design - Chapter IV 

The following examples are given as part of this intro- 
duction to show that separation of two signals with the 
same power spectra is plausible. 

The objective is to estimate the strength of a signal 
in the presence of noise when the noise and the signal have 
identical or almost identical power spectra. In those cases 
(where we can't effectively separate the signals in the 
frequency domain) we will consider the processes of iden- 
tifying the signals in the probability or statistical 
domain. Consider the following: 

Example 1-1 

It is desired to separate a signal from noise when 
they both have a uniform probability density function 
(p.d.f.) of different widths. The amplitude of the noise 
and the signal are unknown and it is assumed that there is 
no correlation between them. The measurement is the sum 
of the signal and the noise. Thus, we assume that: 

1) The noise and the signal are statistically independent 

2) Both have uniform distributions. 

3) The measured signal is given by: 

Z = S + N 
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It is desired to estimate the amplitude of the signal. 
From statistics we know that Z has a p.d.f. which is given 
by the convolution between the two p.d.f., (see Fig. 1-1). 

By inspection of this figure we can see that the p.d.f. of 
Z contains information about the amplitude of the two sig- 
nals. Note that without prior information about or we 
can only determine the two amplitudes, without being able 
to state which is signal and which is noise. The separa- 
tion technique would involve calculation of the p.d.f. of 
Z. Also, this filter can only give us the amplitude of 
the signal and not the whole signal. 

Example 1-2 

Consider the case where the signal is a square wave 
between levels, (Aj^ + A^) and A 2 / where the switching 
time, T, is a random variable with given statistics. The 
noise is the same type of signal with levels + B^) and 
B 2 » Assume that the amplitudes ratio is given: 



A1/A2 



= K, 



B1/B2 



= K, 



( 1 - 1 ) 



Figure 1-2 identifies the signal and the noise. Figure 1-3 
identifies the p.d.f. of the combined signal plus noise. 
With the p.d.f. for Z = S + N and the given ratios and 
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Figure 1-1. p.d.f. of 
their sum. 



two signals and the p.d.f. of 
(Refer to example 1-1). 
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K 2 » one can compute the amplitudes and even determine the 
signal S uniquely as follows; 

Let us assiome that then the estimated values 

of S denoted by S will be (refer to fig. 1-3) : 




S = < 

Ai + Aj 



if or Z^ 

if or 



And for the case 



S 






= < 



A1 + A2 



if Zj^ or Z^ 

if Z 2 or Z^ 



( 1 - 2 ) 



(1-3) 



Notes ; 

(i) In this example there is no assumption one way or 
the other about correlation between the signal 
and the noise. 

(ii) The two cases (Aj^ > B^ or A^ < B^) can be easily 

determined by checking the p.d.f., (the four levels) 
and comparing and K 2 * 

(iii) The levels can be computed easily by the four 
equations ; 



= A, 



B, 
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if 



^2 ®2 ^ 






^2 = 



Ai + A 2 + B 2 if A^_ < 



(1-4) 

if ■^l ^ 

if Aj^ < 

Z4 = Aj + Aj^ + + Bj 

And the cases A^^ > or A^ < B^^ can be determined 
by the two ratios : 




Aj/A^ = 

B1/B2 = 

(iv) The signal can be determined even if the noise is 
the Slim of a square wave plus random noise (for 
example, white Gaussian) . In this case we have to 
introduce the statistics of the switching time 
X (t) and the random noise. 

In the noise free case we can show the processor by a 
simple block diagram given in Fig. 1-4. 
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processor given for exampl 



In these two simple examples , it has been shown that 
in some cases signals can be separated from the noise, or 
the parameters of the signal determined, even if they have 
the same power density spectrum. Thus, the idea of filtering 
a signal in the probability domain seems plausible. 

A second possibility of separation to somehow change the 
measurement, seems reasonable. An example would be in a 
control system where we design our measurement such that 
the system will be observable (see next example) . 

Example 1-3 

Consider the case where the measurement and the state 
equation is given by: 

X(K+1) = X(K) + U(K) 

Z (K) = CX(K) + V(K) 

In the case where the system is not observable there is no 
way to estimate X unless we change something in the system. 
One way to do it is to change the matrix C, so that the 
system will be observable. Let's now consider a practical 
example related to the dissertation. We will consider the 
conical scan radar (see next chapter for details). In c.s. 
we often have the problem of cross talk between the two 
channels (Az and El) • A simplified block diagram of the 
radar is given in fig. 1-5. The system can be described by: 
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Fig. 1-5. Simplified block diagram of 
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e{K) = X(K) - U(K) 

Z (K) = C(e(K) + n (K) ) 



where all the vectors are two dimensional and the matrix C 
is given by: 




e(k) contains the errors in both hn and E^- a and b are 
the cross table factors in the measurem.ent . 

Obviously where a = b = 1 the system is not observable and 
we can't estimate the error (the target position). In the 
case of an observable system, but with a and b close to one, the 

regular Kalman filter estimator converges slowly. In order 
to improve the time constant of the estimator one suggestion 
is to change Z to Z ' by (see fig. 1-5 suggestion #1) : 



Z ' (K) = E • Z (K) 



where 



1 -a 



E = 



-b 1 
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This will lead to a new measured vector (Z'(K)) given 
by; 



Z' (K) = (l-ab) (e(K) + n(K) ) 



This suggestion will solve the cross talk problem but 
it will increase the effects of the noise n . The problem 
will be solved only if a*b 7 ^ 1. So another way must be 
found. This can be done quite easily by changing the modulator 
reference signals. The solution will be that C I and the 
statistics of the new measured noise will be almost the same 
as the given noise and this is indeed an improvment of the 
system performances (see fig. 1- 5 suggestion #2) . The 
details are beyond the scope of this chapter. 

The situation of having the signal and noise with 
approximately the same power spectra occurs in radar. For 
example, we may have; 

(i) Target plus clutter. 

(ii) Target plus target in the same resolution cell. 

The work which follows will concentrate on the tracking 
radar case. We will treat the problem with two major 
approaches ; 

(i) Assume that the radar receiver is given and we have 
to process the output signal. This will give us the 
possibility to try to work the problem using the 
first proposed technique. 
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(ii) Assume that we can change the receiver and also 

make some changes in the antenna such that we will 
get a new measured signal from the receiver which 
will allow us to separate the two signals coming 
from two different sources (targets) . This is the 
second proposed technique. 

Because the second approach will give us the most power- 
ful results, we will concentrate mainly on this approach, 
but we will include all of the information necessary for the 
first approach. 

B. APPLICATION TO TRACKING RADAR 

In radar tracking we can have the real target and the 
false target (chaff, clutter) which are in the same resolution 
cell. (AR<j ^t, Aa<02^j^, where At = pulse length, ~ k»s3m- 

width of the antenna.) Because of the false target, we 
cannot measure the position of the real target 
(angles and range) . For tracking, we must be able to 

measure the error in range, azimuth and elevation where the 
error equals the deviation of the actual target position 
from the estimated position as given by the midpoint of the 
resolution cell. In order to accomplish this, the radar has 
two (at least) windows for each coordinate. The "weight" of 
the target in each window is measured by the difference 
between the normalized weights of each window for each 
coordinate, which is used in calculation to determine 
target position. 
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Example 1-4 



In range we have two gates, the early and late gates. 
By measuring the weight of each gate we compute the error 
in range as : 



AR = R 



^EG " ^LG 
° ^EG ^ ^LG 



R 



= c 



Atr/2 



(1-5) 



where : 



(i) 


AR - 


the computed 


error in range. 


(ii) 


2A'C - 


the width of 


the range gate . 


(iii) 


^o ■ 


the position 


of the range gate. 


(iv) 


C 


velocity of , 


light. 



(V) V 



EG 



/ V(t) dt 



t -At 
o 



t +A-t 
o 



(Vi) V 



LG 



/ 



V(t) dt 



(vii) V(t) - the received signal. 



The time representation of the range gates , received 
signal (video), V_^ and V^_, is given in Fig. 1-5. 

' bo LCjf 
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The same type of equation holds for the angle. In 
this case we have two or more antennas for azimuth, and 
two antennas for elevation. In this technique we estimate 
the weight of the real target in each window, and compute 
the error to close the tracking loop and try to reduce the 
errors to zero. When we have two targets in the same reso 
lution cell, the computed error would depend upon the 
"center of gravity" (radar center) of the radar return of 
the two targets. Because the two targets have the same 
(or almost the same) power density spectra, we cannot dis- 
tinguish between them by spectral analysis; thus we have 
the problem discussed in Section I. A. 

The object of this research is to estimate the true 
position of the target when the false target is near the 
true target. In the research /we apply the following 
practical constraints: 

(i) An estimator which can be used on small missiles. 
We need simple algorithms so that the memory 
size and the amount of computation will be 
reasonable . 

(ii) The size of the antenna is fixed, and we are 

restricted to an antenna of reasonable structure. 
For example we can use a monopulse antenna with 
changes (referred to regular monopulse) but we 
cannot use an antenna with more than four outputs 
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(iii) We will not search for a solution in the r.f. 

range because for each r.f. it is possible to 
have counter measure chaff adapted to that fre- 
quency or frequencies. 

(iv) We will restrict ourselves to surface targets 

which imply no reasonable Doppler shift between 
the two targets. In the problem of air targets 
the difference in velocity between the target and 
the chaff is distinguishable and they can be 
separated easily by a moving target indicator 
(MTI) . 
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II. BACKGROUND THEORY 



A. INTRODUCTION TO TRACKING RADAR 

Before we establish the problem and our solution 
we present a short summary of pertinent tracking radar 
theory. 

The function of tracking radar is to select a particular 
target and follow its course in range, angle and sometimes 
frequency (a direct measurement of velocity) coordinates. 

B. TRACKING BY RADAR 

Tracking systems can be achieved by two different 
techniques . 

Track-while-scan (TWS) 

This method of tracking basically uses data from search 
radar. The idea is to take a sample measurement of the 
target each time that the antenna is pointed to the tracked 
target. This tracking class is not of interest in this 
study. 

Continuous Tracking Radar 

In this method of tracking, the antenna is always 
directed at the target due to function of a control system. 
There are several methods of tracking (i.e., method of 
generating the error signals to close the antenna position 
and the range loops). The important methods are: 

(i) Conical scan (C.S.) 

(ii) Lobe switching (L.S.) 

(iii) Monopulse. 

28 



Those methods are distinguished by their difference in 
angle tracking. The principle of range tracking is almost 
the same for all the radars. Since the most accurate and 
interesting case for us is the monopulse, we will concen- 
trate on it. However the theory of the other techniques 
is summarized. 

1. Angle Tracking 

a. Conical Scan (C.S.). In C.S. systems, angles 
are measured by a single antenna, whose radiation pattern 
rotates periodically about a certain axis. When the tar- 
get is in the axis direction, the radar will have a constant 
return signal. The coordinates of the off axis target is 
determined on the basis of comparison to the envelope of 
received signals with a reference signal. For this reason 
we can not measure the error in one pulse but we have to 
wait at least one period of antenna scan to get the informa- 
tion. A simplified block diagram of the C.S. system is 
given in Fig. 2-2.1. The antenna A, scans in space with 
angular frequency, Q. Two references for the phase detec- 
tor are generated by the antenna (Az . ref. and EL. ref.). 

The position of the antenna is controlled by the control 
system. Assume a target at position T. Because the target 
is assumed to be off the axis of rotation, the echo signal 
will be modulated at frequency Q (the C.S. frequency). The 
amplitude of the signal depends on the position of the 
target with respect to axis of rotation, and the phase will 
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depend on the direction of the angle between the target 
and the rotation axis. When the antenna is on the target 
we will get zero modulation at the receiving signal (see 
Fig. 2-2.2). The modulation pulses pass through a receiver 
which is controlled by an AGC (automatic gain control) . 

We will see later on why the AGC is needed. The signal 
then passes through a "box-car" circuit which changes the 
amplitude modulated pulses into a signal close to sinusoid. 
This signal passes through two-phase detectors (for Az . and 
EL. ) . The outputs from the phase detectors are the errors 
in EL and Az, and by using these errors we can close the 
loop on a control system which will move the antenna to 
reduce the error to zero, so that the antenna "looks" towards 
the target. We can write the amplitude of the signal at 
the output of the receiver as : 



S’ (t,e^) 



S -Gav (6_) [l+m(9_) cos (fit+4) )] (2.2-1) 

o T T c 



where 



S'(t,9^) = envelope of the amplitudes of a 

sequence of radar return pulses. 



S 



o 



the amplitude of the signal (dependent 



on the RCS (Radar Cross Section) of 



the target, the range, the gain of 



the receiver, etc.) 
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1 



the average gain of the antenna, in 
the direction of the target, which 
can be determined approximately by 



Gav(0^) 



G(9^- e^) + G(e^ + e^) 



2 



( 2 . 2 - 2 ) 



to a first order of approximation, 
this average gain is constant, i.e., 
independent of 6^. 

The position angle of the target. 

Conical scan frequency of the antenna. 

Phase of signal relative to C.S. 
frequency. 

The modulation index which is approxi- 
mately equal to: 




S 



max 




(2.2-3) 



2 • S 



av 




where (see Fig. 2.2-2) 



S 



max 



is the maximum received signal 



s . 

min 



is the minimum received signal 




S 



s + s 

max min 



av 



2 



Because we would like to have the same amount of error for 
the same error angle, independent of the target and range, 
we have to measure the modulation index. The method which 
is usually used to achieve this objective is the AGC (auto- 
matic gain control) . The purpose of the AGC is to fix S 

3. V 

such that will be a constant. This can be accomplished 
by controlling the gain of the receiver. The signal can be 
represented by: 



where C is a known constant, independent of the target RCS 
and other parameters. 

If we pass the signal through a highpass 
filter, we will get (after the Box-Car) : 



portional to m(9 ), and phase equal to <j> . With both of these 

X C 

we can determine the error signal after the phase detector: 



S' (t,9^) 



C[1 + m(9^) cos (fit+(})^) ] 



(2.2-4) 




Cm (9^) cos (fit+4)^) 



(2.2-5) 



Now we have a signal with its amplitude pro- 
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J 



AAz 



K*m( 9^) cos(|)^ 



( 2 . 2 - 6 ) 



AEL 



K*m( 6 ^) sincj)^ 



Because m(9^) depends on the antenna parameters, 
the error measured by the system will not be a linear func- 
tion of 9^. The error signals can be determined linearly 
only for small displacem.ent angles. 



achieved by switching the antenna beam between four posi- 
tions (two positions per direction) . We will describe the 
method in one dimension since the two dimensions case is a 
straightforward extension of the one-dimensional situation. 
Refer to Fig. 2.2-3. Because we have two antennas and we 
switch between them, the receiver output belongs to antenna 
1 (S^) and to antenna 2 (S 2 ) • As in the C.S. system because 
of AGC, equal a constant and the error can be 

determined by the difference between and i.e.. 



b. Lobe Switching (L.S.) The method of L.S. is 



error (9^) 



K(S^(9^) - S2(3^)] 



(2.2-7) 



If we have four antennas, by AGC, we keep 



+ S 2 + equal to a constant, so then: 



AAZ 



K[ (S^ + S 2 ) - (S 3 + S^) ] 



( 2 . 2 - 8 ) 



AEL 



K[ (S^ + S 3 ) - (S 2 + S^) ] 
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Figure 2.2-3. L.S. antenna pattern and output signal (one-dimensional only). 
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c. Monopulse Method. In the two previous methods 
(C.S. and L.S.), the measurement of angular error in two 
coordinates (AZ and EL) requires that a minimum of three 
pulses be processed. In practice, the minimum number of 
pulses in L.S. usually requires four pulses (one for each 
antenna). C.S. usually requires much more than four (20 
pulses is a "good" niimber - which is dependent on the fre- 
quency of the C.S. and the P.R.F.). The requirement for 
both of them is that during the scanning or switching the 
signal amplitude does not contain additional modulation 
components other than the modulation given by the antenna. 

If the amplitude contains additional modulation components, 
the tracking accuracy might be degraded, especially if the 
frequency components of the fluctuations are close to the 
C.S. frequency or the L.S. rate. These functions* can be 
caused, for example, by a fluctuating target cross-section. 

In order to avoid this problem, we prefer to base our 
measurements on one pulse rather than many. When a signal 
is received by two different antennas, the difference between 
the two signals might be one of phase and/or amplitude. 

This is the premise of the monopulse radar, i.e., use the 
difference between signals arriving at the same time by the 
antenna. The name monopulse is used to describe those tech- 
niques which derive angle-error information on the basis of 
a single pulse. (Note: In some cases monopulse is called 

"simultaneous lobing.") 
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Two major techniues are used in monopulse: 

(i) Measuring the relative phase of the carrier of 
the I.F. between the pulses coming from each 
antenna. 

(ii) Measuring the relative amplitude of the pulse 
received in each beam. 

Fig. 2.2-4 presents a block diagram of general monopulse 
(i.e., amplitude - phase type), for one coordinate. We 
will see later that by zeroing the distances between the 
antennas, i.e., d = 0, the system reduces to amplitude 
comparison monopulse, and by zeroing the tilted angle, 
i.e., 0^ = 0, it reduces to phase-comparison monopulse. 

The extension into two dimensions is straightforward. In 
most of the cases the sum and the difference are taken with 

4 

the r.f. signal (right, after the antenna), after which 
there are two receivers (one for sum channel and one for 
the difference channel) . This is because of practical 
reasons ; 

(i) Only the r.f. channel receiver must be 
identical . 

(ii) For two-dimensions we save one receiver (3 versus 4) . 

In the analysis which follows we will show the 
sum and difference taken in the i.f. This change does not 
make any difference for monopulse radar but it is necessary 
for the complex angle (CA) method. The signals received 
from the target at antennas and A^ are fed to the receiver 
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Figure 2.2-4. Block diagram of a general monopulse type (mixed type) for one 

coord inate . 



and then the sum (T) and the difference (A) are forming. 

Thus we can write the following cos term of u : 

IF 




5^ + 82 = (0^) cos (ojjpt + (j)) 

+ G2 Ot) cos (ojipt - 4>) ] 



( 2 . 2 - 9 ) 



A( 9 ^,t) = S^-S2 



^A^o ( 9 ^)cos (ojjpt + (}) - 4;^) 



- G2 ( 9 ^) cos (ojjpt - (p - ^^) ] 



where : 



9 

9 



T 



o 



K 



I' 



K, 



Target's position. 

Fixed antenna offset angle. 

Gains of sum and difference channels 
respectively . 



= The average constant signal in the sum 
channel determined by the AGC. 



= The phase shift between the sum channel 
and the difference channel caused by 
the circuits . 



(j) = The phase shift in antenna and A2 

caused by the displacement of target by 
angle 9^ relative to the bore sight 
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direction (see Fig. 2.2-5), and is 
given by: 



2 d 



3^,82 = the I.F. signals within a 

single radar return. 



R + d sin 9^ 

R - d sin 6 _ 
T 



^(R,-R2> 



47Td 

A 



sin 



( 2 . 2 - 10 ) 



Distance between phase centers of 
antenna and (see Fig. 2.2-5) . 



0^,82 = Radiation patterns of antenna and 

A^ respectively; for most of the cases 
they have the same shape only shifted, 
so we can write; 



G^(9) = G( 6 ^ - 9) 
G 2 (e) = G(9 q + 9) 



( 2 . 2 - 11 ) 



where G is a standardized radiation pattern, and 9 is an 
arbitrary direction. 

Let us assume that the gains in the sum and 
the difference are the same (not a necessary condition) : 



K 



I 




( 2 . 2 - 12 ) 
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Figure 2.2-5. Wavefront relationships. 
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By means of the AGC we fix the average sum channel to 
remain constant, as in C.S., we can assume that after the 
AGC we have; 



K • = constant (2.2-13) 

Thus, for any level of signal we will have the same output 
error characteristics. After the phase detector we have 
an output which is the average of the multiple signal 
(^ • A) . It is shown later on that this quantity is pro- 
portional to the error between the direction of the target 
relative to the antenna. 

error (0^) = K' • Average {[ ( 9^, t) * A ( 0^, t) } 

= a'[G^(9^-9^) - G^(9^+9^)]coSii;^ (2.2-14) 

+ [2aG(9^-9^) •G(9^+9^)sin2(J)]sini|)^ 



where; 



K' is a constant determined by the AGC. 
a' is a known constant. 

From the last equation we can see that the 
error is only a function of 9^ and equals the deviation 
from the bore sight (the rest are known constants determined 
only by the system and not by the target) . So from 
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we can choose the parameters d and 6^. In the two major 

types of monopulses, we choose these parameters: 

(i) Amplitude comparison monopulse: 

Choose = 0 or d = 0 . Practically, we choose both 
of them equal to zero. We would get, for the case 
= 0 and d = 0 : 



(ii) Phase comparison monopulse: 

Choose = tt/2, or 9^ = 0*. Practically, we choose 

both. Then, for the case 6 =0 and ^ = tt/2: 

' o o 



error (6^) 



[G^ (9 -e^) - (9 +0„) ] 

t. ' o T o T 



(2.2-15) 



error (9^) 



2aG^(9^) sin 



2aG^(9^) sin(ysin9^) 



(2.2-16) 



Y - 8TTd/A 



From the last equation we see that 9^ can be 
determined. Near the origin (®rp ^ 0)* Then we 



can assum.e: 



Assuming G(0) = G(-0). 



44 



(2.2-17) 



G(0 



T' 



G 

o 



and 



sin(T ‘Sin 6^) :: Y 9^ 



so that: 



0^ = 3-error(0^) (2.2-18) 

where 3 = known constant, where the error is the 
output after the phase detector 
(3 = 1/uG^^Y). 

When we wish to determine the position in two 
coordinates, theoretically we need only the addition of one 
antenna; practically, four antennas are used. The theory 
is almost the same; the only difference is that now we 
have three channels: 



(i) 


Sum channel 


_ V 






(ii) 


Difference 


channel 


in AZ 


\z* 


(iii) 


Difference 


channel 


in EL = 


'^EL' 



where (refer to Fig. 2.2-6): 

I = + S2 + 

= (S^ + S^) “ (S 2 + S^) (2.2-19) 
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Figure 2.2-6. Relative location of the antennas 

in two-dimensional monopulse. 
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and 



= the output from antenna i. 

2 . Range Tracking 

The most widely used technique for tracking in range 
is based on two range gates. There are no special tech- 
niques as for the method of angle tracking. Information about 
the range is taken in most of the cases from the sum 
channel in the monopulse radar and from the incoming 
pulses in C.S. or L.S. We will briefly cover the prin- 
ciple of automatic range tracking. A block diagram for 
the range loop tracker is given by Fig. 2.2-7 and the 
time representations of the signals are given in Fig. 2.2-8. 
We have two major signals in the loop: 

(i) The video pulses - S(t) 

They contain information about the range position 
of the target. 

(ii) The range estimate - R. 

The main idea is to open two gates, g^/ before and 

after the estimated time of the center of arrival of the 

2 . R 

pulse (t = , R = taraet range, c = velocity of the light) . 

c 

One is the early gate (g^^) , and the other is the last 2 ; gate 
. The portion of the signal contained in the early gate 
is subtracted from the portion of the signal contained in 
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Figure 2.2-7. Block diagram of tracking in range. 
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Figure 2.2-8. Time representation of the signals in 

Fig . 2.2-7. 
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the late gate after integration. The error is calculated 
after the end of the late gate. The magnitude of the error 
signal (^R) is a measure of the difference between the cen- 
ter of gravity of the video signal, S(t), and the estimated 
range, which is the center of the two gates. The error is ' 
fed into a control system to estimate the range. The output 
of the control loop is the estimated range, R. To convert 
the range into time (to generate the gate's pulse) we use a 
converter which has the time of the transmitting pulse as 
a refernce input. The outputs of the converter are: 

(i) The two range gates (g^ and 
(ii) A reset pulse to reset the integrator before the 
range gates. 

Because we want to have the same output error, Ar 
for all the targets, independent of the target, we must 
normalize the pulses, and this is done by AGC, which has 
already been applied in angle tracking. 
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III. STATEI-IENT OF THE PROBLEM 



A. INTRODUCTION. INTERFERENCE CAUSED BY TWO TARGETS 

Let's consider the problem of tracking a target when 



a false target is 


in the same resolution cell, i.e.. 


AL 

R 


< 9 

3db 


AR 


< At- c/2 



where : 



®3db 


the beamwidth of the antenna (in radians) 



2At = pulse length 



R 


range of the target 


AL 


the difference between the true target 
and the false target perpendicular to 
the line of sight of the antenna 


AR 


the difference between the true target 
and the false target in range (parallel 
to the line of sight) 


c = 


the velocity of light. 



For simplicity, let's consider the two-dimensional 
case only (range and elevation) . The principle in three 
dimensions is the same but more complicated, and the results 
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can be obtained by a straightforward extrapolation of the 
two-dimensional case. We will consider the case of point 
targets for developing the theory, but we will implement 
the results when the targets have complex structure. 



B. INTERFERENCE IN ANGLE CREATED FROM TWO TARGETS 

Consider the case of two point targets (or spheres) 
in the same resolution cell of the radar with spacing 
AL (i.e., ^ ^ velocity of light). 

As the relative path lengths between the radar antenna 
and the two sources vary (i.e., R^^ and R 2 ) / the two sig- 
nals will alternately add and subtract, and so the amplitude 
and the phase of arrived signal will change. Although such 
a simple situation (two points targets) may be fictitious, 
it will illustrate the main behavior. The relative ampli- 
tude between the RCS of the two targets is assumed to be a 
constant "a" and the relative phase difference also constant, 
a. The difference in phase is due to difference in range 
(Ar) or to reflecting properties. The relative angular 
error (A0/6p) is given by [23]: 




+ a cos a 



1 + a + 2a cos a 



( 3 . 2 - 1 ) 



The position of the stronger target corresponds to 
A0 = 0, while the smaller target position is at A0/0^ = 1. 

(See Fig. 3.2-1.) The position of the tracking system depends 
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can be obtained by a straightforward extrapolation of the 
two-dimensional case. We will consider the case of point 
targets for developing the theory, but we will implement 
the results when the targets have complex structure. 

B. INTERFERENCE IN ANGLE CREATED FROM TWO TARGETS 

Consider the case of two point targets (or spheres) 
in the same resolution cell of the radar with spacing 
AL (i.e. , ^ < 63 ^/ ^ velocity of light) . 

As the relative path lengths between the radar antenna 
and the two sources vary (i.e., R^ and R 2 ) , the two sig- 
nals will alternately add and subtract, and so the amplitude 
and the phase of arrived signal will change. Although such 
a simple situation (two points targets) may be fictitious, 
it will illustrate the main behavior. The relative ampli- 
tude between the RCS of the two targets is assumed to be a 
constant "a" and the relative phase difference also constant, 
a. The difference in phase is due to difference in range 
(Ar) or to reflecting properties. The relative angular 
error (A0/9j^) is given by [23]: 



2 

A0 _ a + a cos a 

9 2 

D 1 + a + 2a cos a 



(3.2-1) 



The position of the stronger target corresponds to 
A0 = 0, while the smaller target position is at A 0 / 0 j^ = 1. 

(See Fig. 3.2-1.) The position of the tracking system depends 
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on the relative phase a and the ratio a . One can show 
that for 0 £ a £ 1 and 0 £ a £ 2tt - the values of A9/A6j^ 
will be between 



— oo < 




0.5 



(3.2-1) 



When the echo signals are in phase (a = 0) , 
the error reduces to — - j^ / which corresponds to the so- 
called "center of gravity" of the two targets. 

Now, when we have a complex target, i.e., the 
ratio a is a random variable which changes from pulse to 
pulse ^ Also a is random variable which changes from pulse 
to pulse _ Taking these into the given servo loop of the 
radar system, one can compute the statistics of the random 
variable A0, and this can give us a good approximation of 
the error (for example, mean and variance) caused by the 
addition of two targets in the same resolution cell. How- 
ever, the simple case results are sufficient for our problem, 
i.e., to show that the result of two target returns in the 
same resolution cell causes an error in the estimation of 
the position of the true target which depends on the signal 
ratio between the two targets . When the two targets have 
a complex structure rather than a point, the results are 
much more complicated and we have to take into account the 
circuits involved and the statistics of the target's returns. 
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C. INTERFERENCE IN RANGE CAUSED BY TWO TARGETS 

In range measurements we have a similar effect, called 
glint, i.e. , we will get an error range depending on the 
difference of range between the targets and the phase 
difference. We can write the measured error in range as: 

t^+Ax 

o o 

AR = K[ / S(t)dt - / S(t)dt] (3,3-1) 

t 

o o 

where 

S(t) = S^(t) + S2(t) 

2 At = the width of the range gate 

t^ = the estimated range of the target (i.e., 

the position of the range gate) . 

The behavior in range is strongly influenced by the 
radar parameters and the circuits used for measurements 
(for example. At), but a "good" estimation of the behavior 
can be taken as when the two targets are totally in the 
range gate (i.e., the two targets are very close in range, 
or there is a large range gate compared to At) , so that the 
range gate will be at the central gravity of the two targets. 
Also we will assume that the video pulses are rectangular. 
Refer to Fig. 3.2-1; the center of gravity of the two 
pulses will be: 

T^g = j t s(t) dt / J s(t) dt 

ft] [t] 
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Figure 3.2-1. Simplified time representation of the two pulses. 
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Amplitudes of the two signals, separately. 

The two signals combined in phase. 

The two signals for combined out of phase. 

Received signals for the general case. 

The center of gravity of the target. 

■ The center of gravity of the two combined pulses. 

Figure 3.2-2. Range errors caused by two targets. 
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s(t) dt] 



(3.3-3) 



, T^+T„-2At 

T = -{ = 

eg 2 



T^-At 



T^-At 



+ [ 



T 1 +T 2 



T^+At 



S(t)dt] + [ 



T^+T2+2At 



T^+At 



T2~At 



S (t) dt] } 



Tj_+At 



where in case of < T^; and 2At - AT > 0; 



S(t) =< 



t < T^-At 



; . T^-At < t < T2~At 



a* ; T 2 ~At < t < T^+At 



T^+At < t < T2+At 



0 ; t > T2+At 



B ^ 



T2+At 



S (t)dt 



T^-At 



I * = 



ai t a2 e 



3« 



a is the relative phase between the two 
returned pulses. 



In general where is the position of the true 

target. Thus we get an error in the computed range. This 
error is between the two extremes given by the cases when 
the two pulses are combined in phse and out of phase. When 
the two pulses are combined in phase we get a simple 
expression for (3.3-3): 
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( 3 . 3 - 4 ) 



T 

cgs 



ai + a2 



The error due to the addition of the second target is 
given by (for the case that they are in phase) : 



AR 



s 




^2 



ai + a2 




( 3 . 3 - 5 ) 



which is not equal to zero when ^ '^2* 

/ 

D. CONCLUSION 

In the last two paragraphs we have shown that when 
there is an additional target in the resolution cell, the 
radar will not track on either of the targets, but on 
their center of gravity. In order to eliminate this 
phenomena, we must change the design of the radar so that 
it will track on the two targets. 

There are two main ways for accomplishing this objec- 
tive which have been developed in this study. 

(i) Extract information of the location of the target 
by processing the signal after the recevier. This 
leads to the implementation of what we have called 
the probability filter (Chapter V) . 

(ii) Modify the receiver and the filtering process after 
the receiver in order to extract the information of 
the two targets. 

This leads to an antenna modification and a new type of 
modified Kalman Filter (Chapter IV) . We will consider later 
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the case of a ship as the true target and chaff as the false 
target. This implies very small Doppler shift between 
the two targets. In the problem of air targets the differ- 
ence in velocity between the target and the chaff is 
distinguishable and they can be separated easily by a 
moving target indicator (MTI) . In summary, the following 
techniques are presented in this dissertation: 

(1) Antenna modification plus a Modified Kalman Filter 
(M.K.F) . In this case the location of the targets 
are resolved on a pulse to pulse basis (see Chapter 
IV-A) . Thus it can include the case of air targets 
where doppler shift is not applicable rather than 
be restricted only on surface targets. 

(2) Probability filter . This solution is suggested in 
case of C.S. or L.S. radars. Here the location of 
the targets cannot be resolved on a pulse to pulse 
basis (see chapter II) , and many pulses must be used 
in order to solve the problem. Because the dynamic 
of a surface target is slow, we can assxime a 
stationary process for the target's RCS. Thus a 
time average over a small interval can used as a 
good estimate for the mean of a function of the 
RCS. 
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V. A MODIFIED ANTENNA AND KALMAN FILTER SOLUTION FOR 
UNRESOLVED TARGETS - FOR A MONOPULSE SYSTEM 



A. THE THEORETICAL SOLUTION 
1 . Introduction 

The problem of two targets in the same resolution 
cell is discussed in the literature, mostly from the point 
of view of resolution and multipath [Refs. 12,13,14,15,16, 
17, 18, 19, 20]. The principles involved in the multipath 
problem are similar to the two-target problem, and it is 
proposed that they be adopted for this case. The major 
difference between the two problems is that in the multi- 
path problem the amplitude ratio and the relative phase 
between the two return signals are known. This information 
is missing in the two-target problem. The method used to 
solve the multipath problem is called the "Complex Angle" 
method (CA) . In this chapter we resolve in a new way the 
positions of the targets in a single pulse and demonstrate 
how to filter the data. We break the problem into two 
parts : 

a. Resolve the positions of the targets and their 
amplitudes, assuming that there is no noise. This 
gives us the positions of the targets plus estimation 
noise . 

b. Filtering the positions of the targets. Since the 
foregoing yields the position of the target plus 
estimation noise, an additional filter must be used. 



I 



i 






Furthermore, the estimate of the positions is with- 
out the decision of target vs chaff. Thus this 
filter must select the two positions related to the 
chaff and the target. 

The general assumption is that in the input to the antenna 
there is a white Gaussian noise. Since the first block 
includes a non-linear transformation of the signal (see 
fig. 4.1-1), there is non-Gaussian noise at the output. 

This output is the input to the filtering block. Our 
filtering procedure is complicated even with the assumption 
of a Gaussian process. Thus, to save a lot of computation, 
the assumption is made (similar to the usual monopulse analysis) , 
that all the variables are Gaussian. In regular monopulse 
systems the measurement process requires multiplication 
(see II. B) . The variables A and I are assumed to be 
Gaussian. Hence the measured error is not a Gaussian varia- 
ble although in most of the regular monopulse radar, is 
assumed to be Gaussian variable. The same type of assumption 
is made in this application. 

The block diagram is given in Fig. 4.1-1. We 
start with the first block, i.e., the resolving procedure. 

2 . Resolving Targets with Monopulse System 

This paragraph introduces a new technique for 
resolving two targets, based on the "complex angle" method. 

The uniqueness of this solution is that we use only a four 
element antenna, which is the regular monopulse. However 
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we require four output channels from the antenna rather 
than three channels as in regular monopulse. The resolved 
positions, similar to regular monopulse, are made in one 
pulse . 

For monopulse system and two targets, the output 
of each antenna is : 



®i ^ ■^l^i ^^I'^l^ ^2*^i ^^2 '^2^ 1,2, 3, 4 (4.1-1) 



or, in vector form; 



where 



e = A^G(X^) + A 2 G(X 2 ) 



G. 

1 



th 

the voltage pattern of the i received 
beam 



X 



-a 



A 



a 



coordinates of target a; (X^,Y^) = 1,2 

effective complex voltage due to target ^ 
(normalized to the range) , and is given by 



^ (4.1-2) 



(ju^ = the doppler frequency of target a 

CKX 

9^ - phase associated with target' a 



a 



th 



a - The a target voltage cross-section 
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I 



- one-way voltage pattern of the transmitted 
beam 

K - normalized factor. 

Note: For implementation, it is much simpler to assume 

= 0 and then solve the equations. This is reasonable 
if we deal with surface targets, and so 0, and this 

assumption would simplify our decision procedure. However, 

the reason that we do not assume = 0 is because 
Yj^ is a random variable with almost zero mean but with 

variance not equal to zero. For example, a ship is tilted 
by the waves and the c.g. of the ship is changing with 

time. Thus, the assumption that Y^ = 0 can cause an 

error. Also, the assumption Y = 0 would lead to a special 

case solution which is not desirable. 

There are four unknowns for each target: 

The positions - X = (X ,Y ) 

^ ~a ' a a 

The amplitudes - A = (Re A , Im A ) 

^ ~a ' a a 

Because eq. 4.1-1 is a complex equation and we have 
four antenna elements, we have a total of eight equations. 
After solving (4.1-1), one can use (4.1-2) to solve for 
and the angle • By using at least two different 

samples, we can also separate from 0^ and have the 

doppler information about the target. The main approach of 
this section is the solution to (4.1-1). Because this 
equation is quite complicated, the main idea is: 
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( 1 ) 



Solve it by approximation, using Taylor series [12] 
as is usually done in regular monopulse analysis. 

(2) Solve the non-linear equation with -a starting point 
given by the first step. 

The Taylor series solution for this case follows 
with the solution given by (4.1-18,19). 

Our assumption is that the beam is composed 
of factors, each related to a particular angular 
dimension, i.e.: 



We solve the equations for the unknowns by expanding the 
gains, and by Taylor series up to first order about 



Gi(X) = X^(x)-Y^(y), i = 1,2,3,4 



(4.1-3) 



the boresight line, i.e.: 



X^(X) 



/g (1 + a^X) + 0 (X^) 



(4.1-4) 



y^(y) 



^ (1 + 6^y + O(Y^) , i = 1,2, 3, 4 



where 



g - is the boresight gain 




is the slope of the beam near the boresight 



line for X direction 



8^ - is the slope of the beam near the boresight 
line for Y direction. 
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From (4.1-3) and (4.1-4) we get: 



G^(X) = g(l+ + 6 j_Y + a^B^XY) 



(4.1-5) 



i = 1 / 2 , 3, 4 



Equation 4. 1-1 can be approximated by: 



= gA^ (1 + a^Xj^ + Bj_Y^ + a^Bj_X^Y^) 



(4.1-6) 



+ gA 2 (1 + a^X^ + B^Y^ + a^Bj_X^Y^) i = 1 , 2 , 3, 4 



Those are four complex N.L. equations. Let us now 
"linearize" the equations by substituting new variables; 



= A^ + A2 



^2 ^ 1^1 ^ 2^2 



(4.1-7) 



Z3 = A^Y^ 4 - A2Y2 



Z 4 = Vl^l + ^2''2''2 



In this case the linear approximation of (4.1-6) is 



e^ = (Z^ + a^Z 2 + 6 j_Z 3 + a^B^Z^) -g , i = 1,2, 3, 4 (4.1-8) 
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or in short form; 



e = B • Z 



(4.1-9) 



where ; 



Z ^ 



- ^4 



B ^ g 



1 



# • 



1 ct. 



“l^l 



6 . ct,i 6 . 

4 4 4 



(4.1-9a) 



We can assume that B ^ exists and so the solution of Z is 
given by: 



Z = B ^ e 



(4.1-10) 



From the knowledge of Z we have to find the location of 
the targets (X), and the amplitudes. By equation 4.1-7 it 
is easy to determine the amplitudes, given the locations: 
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(4.1-11) 



Z2Y2 - X 2 Z 3 


Z1X2Y2 - 


X1Y2 - X^Y^ 


X2Y2 - X^Y^ 




-^iVl ^ ^4 


X1Y2 - X^Y^ 


X2Y2 - X^Y^ 


1 new variables 


T A 

(a = [a^,a 2 ,a 3 ,a^] ) 


1 — 1 

1 

CN 

CN 

X 
1 — 1 
X 




X2Y2 - X^Y^ 




X2 - x^ 




X2Y2 - 




Y^Y 2 (Xj - X^) 




X2Y2 - X^Y^ 




^2 - ^1 





(4.1-12) 



X2Y2 - 



This gives us a new linear form of equations 4.1-11: 



(4.1-13) 

The last two complex equations (4.1-13) are four equations, 
and we have to write Z,- as it is the sum of real and 
imaginary parts . 



Z^ai + 

^ 1^3 ^ 4^4 
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= Ri + j li . 



i = 1/ 



(4.1-14) 



Let's define 



R 



A 




I 





where ; 



Z = R + j I 



We now can return to eq. 4.1-13 by defining: 



H, 



0 



R^ 



o 
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Then we write equation 4.2-13 in short form; 



H. 






(4.1-15) 



Solving this equation gives 






R1I4 - 



«2l4 - ''4^2> 

(R3I4 - R4I3) 

(-R3I1 + R3I3) 



(4.1-16) 



-1 



Assuming exists, then a is determined. Now we have to 



solve X^. First we notice that; 



X 2 - 



Y - Y 
2 ^1 



X 1 X 4 



(4.1-17) 



Y Y 
14 



Xo - X3 



X 2 Y 2 



XiYi 
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For convenience we define a new variable: 



U = H X 
u~ 



X ^ 



X , 



(4.1-18) 



H 



u 



A 



1 

1 


1 o ' 

1 


0 


D 


-1 


1 

1 — 1 


1 — 1 

1 


1 

1 — 1 



The inverse transform is given by: 



X = ^ H U 

2 u ~ 



(4. 1-19) 



Solving (4.1-17,18)- for U gives, by eliminating x: 
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a 



2 



U 



3 



a 



4 




U 



(4.1-20) 



S 




a 



a 



2 



3 



* - ^ 1^4 + ^ 3^2 



when S is a variable which can be +1 or -1. In determining 
U one solves for U 2 and first. Using these values one 
can find and This is illustrated in Fig. 4.1-2. 

The variable S occurs because we have a nonlinear 
equation. It is easy to see that changing S from +1 to 
-1 will cause the whole solution to change indices. (I.e. 
given e 4.1-10 yields Z, R, I and b 2 / 4.1-16 yields a, 
4.1-10 gives U, 4.1-19 gives X, 4.1-11 gives and A 2 ; 
see figure 4.1-2): 

When S = 1 



X 




A 



1 



A 



A 



2 
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Antennas Synchronized 







Figure 4 . 1.- 2 . Block diagram of determining the 

amplitudes and the positions of 
the targets. 



or when S 



-1 




From which it is clear that in both cases we have the solu- 
tion of the location of the two targets. For the solution 
(n ) assume an arbitrary that S = +1. 



solution 
for S=1 






Then we have a unique solution for the positions of the 
targets but we can't match the position of the target to 
the measurement. There are two possibilities for target #1 
position (and the same for target #2), i.e.: the position 

X X 

of target #1 can be ( or ( ^) . 



( 

i 

I 

I 



The block diagram of solving the positions and the 
amplitudes of the two targets is given in Fig. 4.1-3. 
Everything is straightforward in the block diagram except 
the amplitude vector for the range, which will be explained 
in the following paragraph. 

Solution for Range 

As was mentioned in Sec. II. B, the error in range 
is given by: 



AR 



K 



V - V, 

R V + V. 

eg ig 



(4.1-21) 



where : 



V - is some measure of the amplitude of the 
xg ^ 

X gate (early or late) 



K„ - constant. 



In our case we cannot take the amplitude in each range 
gate and substitute into the formula because at each gate 
we have the sum of two amplitudes. Thus we have to use the 
technique developed for resolving the angle. The simplest 
way is to use the procedure for solving the amplitudes for 
each target, but now we have to consider two amplitudes for 
each target (for the early and late gates) . The procedure 
is described in the block diagram 4.1-3. 

(1) For solving the regular araplitude we use 
formula 4.1-11. We need the delay because we need to wait 
for the solution for X and Y. 
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Figure 4.1-3. Block diagram of the amplitudes. 



(2) We take (at least) two samples of 7. in the 
late gate and in the early gate (Zj^/Z^) and we use the 
same formula for obtaining the amplitudes at each gate, 
where ; 




For the purpose of tracking in range we do not need the 
complex signals. Thus we give the output only for the 
absolute value of each amplitude. 

Note that here we cannot use integration as sug- 
gested in Sec. II. B because we have to take a sample of 
the amplitude. To perform a process like integration we 
have to take many samples . Our process involves a phase 
detector plus a computation for each sample, and this limits 
the samples that we can take. In our case in order to save 
the amount of computations, only two samples are taken. 
Furthermore, we can save even more calculations by assuming 
the range of the target and the chaff to be the same. 

After we have the approximate solution for the 
unknowns, one can substitute these values in eq. 4.1-1 as 
a first approximation and iterate the N.L. equations for 
a more accurate solution where we can assume that: 
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\ 






U^(u) 



K. 

1 



sin C.(u - 6^^) 

q(u - a^.) 



i = 1, . . . , 4 (4.1-22) 



with: 

U. - is X. or Y. 

1 IX 

u - is the coordinate x or y 

- parameters of antenna i 
9 . - the offset of antenna i 



Let us now observe that this method gives us exactly the 
same solution for the single target case. In this case, 
A 2 = 0 and so (4.1-7) reduced to: 



And so 



Z, = A, 



AiXi 



AiYi 



A^X^Yi 




(4.1-23) 



(4.1-24) 
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Now, let's assume the same beam shape for all the antennas 
(i.e., |a^| = \Q^\ = a) . Then the matrix B becomes 

(refer to Fig. 2.2-6 and equation (4.1-9a): 



B = g • 



1 

1 

1 

1 



a 



-a 



-a 



a 



-a -a 



a 



a 



-a 



-a 



a 



a 



(4.1-25) 



Solving for B 



a 



a 



a 



-a 



B 



-1 



4ga^ 



a 



-a 



a 



a* 



-a 



-a 



-a 



a 



(4.1-26) 



-1 



-1 



Now we can solve for (4.1-19): 

2 ^ (e^ t ^^) ” ^^1 ^ ^3^ 
1 a e^ + S 2 + 




a T 



^1 



1 (ei + e^) - (e3 + e^) 



1 

a 




(4.1-27) 



which is the well known solution for regular monopulse 
systems. Thus, from this result we can say that regular 
monopulse system is a special case of our development. 
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3 . The Filtering Procedure 



In the last paragraph we investigated a way to 
resolve the unresolved targets where we have a monopulse 
system with phase detector. The development was based on 
the assximption that we have no noise. We could work this 
problem as was suggested by [20] but this leads to a lot 
of computations, and the estimator functions on a pulse to 
pulse basis ([20] does not use the prior information). Thus, 
after our estimator we have to use a filter. Furthermore, 
we have also to make a decision about target vs chaff. 

For simplicity we assume that the outputs n (see 
Fig. 4.1-1) are Gaussian variables, so we must solve two 
problems : 

(1) Selection - i.e., at each pulse we have to 
distinguish the origin of the measured data, which means 
that each measurement gives us two sets of data, target 
and chaff, but we cannot distinguish between the two. 

(2) Filtering - After we make the decision, our 
data is not noise free. Furthermore, the decision may 
be wrong. For this purpose, we have to filter the data 
output from the selection block. In the previous analysis 
we assumed that measurement noise does not exist. In 
reality we should have started with the equations: 

®i “ Aj^G^(Xj^) +A2 G^(x 2) + i = 1,2, 3, 4 (4.1-28) 
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e 



I 

I 

I 

1 

( 



where is a complex white Gaussian noise: 



^i = 



n . + j n . 

Cl SI 



(4.1-29) 



n . ,n . - white Gaussian noise with zero 

Cl SI 

mean, each of them with variance 
2 

a , assuming no correlation 
between n and n, . 



However v/e used (4.1-1) so that after our processing, we 
can assume that a Gaussian noise is introduced. The 
computed information needed for Gaussian variances are: 

i) The means 

ii) The variances 

iii) The correlation between the noises. 

That is, we have to include, for the random process n: 



y = E{n} 



(4.1-30) 



= COVAR{n} (4.1-31) 

which is straightforward. 

Once we have n output of the resolved block positions with 
the approximated statistics of the noise, we can make our 
decision and filter using the information about the ampli- 
tudes of the signals and the dynamics of the targets. 

The decision and filtering process is described in the 
block diagram given in Fig. (4.1-4). 
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Figure 4,1-4. Block diagram of the decision of target- 

chaff and filtering. 



Our information comes from the monopulse antennas, 
and we compute the positions and the amplitudes of the 
targets in the resolving block. The outputs from the 
resolving block are n and which include the posi- 

tions and the amplitudes of targets number one and two, 
and we have for the filtering procedure to decide if 
belongs to the target or the chaff. The information needed 
to decide is the predicted position of the targets and 
the covariance matrix of the error of our prediction. 

After the selection block, we get the information 
needed for the filtering, i.e., the amplitudes and the 
positions of the target and the chaff. The output of the 
filtering block is the estimated position of the target 
and may be, if needed, the position of the chaff and the 
velocities . 

There is also an outside output to the filtering 
block, i.e., during the process we change the position of 
the antenna in order to zero in on the target. Because our 
model refers to the absolute coordinate we have to take 
the change of the relative position antenna-target into 
account. In our case we make the above assumption for 
simplicity: 

i) The process is sequential. This means that we 

calculate the estimator output and the covariance 
matrix from pulse to pulse. We do not have to store 
all the measured information so far. This assump- 
tion comes from the desire for a simple procedure 
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and reduces the amount of storage and computation 
needed for a non-sequential process, 
ii) Our variables are Gaussian. 

The theory of the selection and filtering is given 
in Appendix A. We give here a short description of the 
procedure. 

i) The dynamic model: 



X(K+1) = 'l>X(K) + AU(K) + W(K) 



X(T) 



(K) 



X(K) ^ 



X 



(C) 



(K) 



= the amplitudes, (4.1-32) 
A and the 
positions of 
the targets , 





1 

o 

1 


o 


1 

1 4(0 

1 







$ 



(i) 



the transition matrix of each target 
(i=T or i=C) 
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ii) The measure model: 



n (K) = L(K) • C • X(K) + V(K) 



C 



[C 



(T) 



C 



(C) 



n 



( 1 ) 



(K) 1 



n (K) 



A 



n 



( 2 ) 



(C) 



L(K) - is the random switching matrix 



(see eq. A-1) . 



iii) The statistical information about the states and the 
noise is given in Appendix A. 

a. The Selection Procedure. All the procedure is 
described in Appendix A. Hence only the block diagram is 
given here (4.1-5). 

Description: We have in this block the above inputs: 

i) The resolved information about the targets - n (K) 

ii) The estimate predicted positions and amplitudes 



of the targets - Y(k|k-1) 

iii) The covariance matrix of the error of the estimate 
prediction - A. 

The output from this block is the assumed information about 
the target. The resolved information of the target is 
transferred into two forms: 
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n^ure 4.1 5. Block diagram of the selection procedure. 



I 



i 











n 



( 1 ) 



i) 




which is the original input. 



n 



( 2 ) 



n 



( 2)1 



ii) 




which is the switched input. 



( 1 ) 



n 



After we form the two possibilities t^, we subtract the 



predicted measure Y(K|K-1) with the formed t^ and we get 
a two error vector r^. We now form the likelihood test, 
and we decide t^^ or t 2 according to > %2 ^2 ^ ^1' 

The output from this block is the assumed information about 
the target. 



made, we have to filter the output from the selection block 
(4.1-5). The filtering could be straightforward if we have 
high confidence that our selection is correct. In our case 
we are not sure if our selection is correct, and we must 
take this into account. Since we cannot tell if our selec- 
tion was correct or not, we can only take into account the 
probability that the selection was correct, i.e. , we cannot 
reject or accept the data, we can only "weight” the data 
depending on the computed probability that our data is 
correct or not. 



b. Filtering. 



After the decision about target vs chaff is 
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The two possibilities are; 



i) A true data available, i.e., we decide the right 
selection. 

ii) A false data available, i.e., we decide the wrong 
selection. 

We must assign a certain "weight" (depending 
on the probability that this data is true) to the data. 
Because all the procedure is descibed in Appendix A, 
only the block diagram is given here, (Fig. 4.1-6). 

B. SIMULATION OF THE MODIFIED KALMAN FILTER 
1. Introduction 

The M.K.F. is very important in our system, and 
we must show that the concept of the modified filter works. 

We will use only one component of measurement to show how 
it performs. Because the most distinguishable component 
between the chaff and the target is in the height, the 
simulation is only for height. To take into account the 
variance and the bias of the estimation error the probability 
of missing is used as a criterion. The probability of 
missing includes both the bias and the variance of the 
estimation error. Furthermore the radar system is assumed 
to be in a missile. 

We have chosen only one component of measuring, 
i,e. , the height. The target which is assumed to be on the 
surf ace, has zero net velocity in this direction. The chaff 
has no negative velocity in height. Thus we can define, 
according to (4.1-32); 
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Figure 4.1-6. The block diagram of the filtering procedure. 



I 



where: 



X(K) = 







X^”^^ (K) 




X^(K) 






1 X^"^^ (K) 


X,(K) 


lll> 


(C) 


<011 


t » 






X. (K) 






(K) 




X 

X2^^^ (K) 




X^^^ (K) 1 

-V. ( 

! 






(4.2- 



X^(K) 



( 'T’i 

(K) 



is the height of the 
target at time K 



X 2 (K) 



(C) 

X^"-^ (K) 



is the height of the 
chaff at time K 



X^(K) = X2^^^(K) 



is the velocity of the 
chaff at time K 



and 

A = 0 



According to the above assumptions, we can relate (4.1-32), 
(A-2) and (A-2a) with the proper matrices- The matrices 
are changed according to the simulations performed. 

Note: In the simulation we use the terms (refer to (4.1-32)) 

"without switching" which means 

L(K) = I for all K's 

"with switching" which means that L(K) is 

as defined in (A- lb) . 
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In the simulations we compare the M.K.F. "with switching" 
to the R.K.F. "without switching" as a reference (ideal) , 
and to the R.K.F. "with switching" to show the improvement 
that we get. 

A Monte Carlo simulation with 200 ensemble members 
is performed. The results are given in comparison to the 
regular Kalman Filter (R.K.F.) (where there is no switching 
in the input). The criterion for missing is as follows; 

The missile is simulated as an inertial mass 
and it is guided to reach the estimated height 
of the target, i.e.; 

Xj^(i+1) = KXj^(i) + (l-K)X^'^^ (i+1) (4.2-2) 

where 

Xj^(i) - the missile's position at time i 

^ (T) 

X (x) - the estimation of the height at 
time i, see (4.2-1) 

K - transition constant of the missile. 

It is chosen to be 0.8. 

The criterion of missing: 

|Xj^(TTG) - X^"^^ (TTG) 1 £ S S ^ 0 



where : 
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TTG - time to go (chosen to be 40 sec) 

(i) - the height of the target at time 

i , see ( 4 . 2-1) . 

In our case 3 is chosen to be 0.75 which experi- 
mentally gives a probability of missing of about 

0.5. 

We study the performance of the two filters by changing 
three significant parameters: 

variance of the noise. 

- mean of the initial states. 

- standard deviation of the initial position of 
the chaff, i.e. : 




[X^CO) - X2(0)]^} 



Note: In the simulation instead of the formula (A-13) for 

a(K) the following formula was used; by mistake: 



where : 



a (K) 



1 

1 + Y (K-1) 3 (K) 



Y(K-l) 



-L/2 

P 



L = b'^(K)A ^(K)b^K) 



T 

b^ (K) 






X;l(k|k-i) 
X^ (k1k-1) 



X2 (k1k-i) 
Xj_(k|k-i) 



and 3(K) is the same as 3(K) given in (A-13a) , and 
A ^(K) is defined in (A-6). This leads to an a(K) 
slightly higher than the real a(K). Apparently this 
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change in a(K) will not affect the simulation results 
significantly. Hence it should not change the 
simulation conclusion. 

2 . Simulation Results 

a. Influence of the noise on the Performance of 
the Filter 

The following parameters and matrices are chosen: 



where 




0 

1 

1 



C = 



1 0 
0 1 



0 

0 



R = 



1 

0 



0 ~ 

1 



M 



50 0 0 

0 5000 0 

0 0 0.5 



X(0) 



10 

100 

-1 



Q 



a 0 0 

0 a/100 0 

0 0 a/4 



"a" - is a parameter to be changed. "a" 

is changed from 0.1 to 100. 
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Figure 4.2-1. Probability of miss as a function of -lo . 



Three results are given here: 

(1) . The probability of missing as a function of -10 lOg (a) 
(Fig. 4.2-1). From the simulation results we see 
that ; 



(i) The R.K.F. (without switching) gives a better 
result than the M.K.F. (without switching), 
but they are close. Thus the MKF approaches 
the ideal. 

(ii) The M.K.F. gives significant improvement over 
the R.K.F. where the input to both is switched 
randomly. (The probability of a miss for the 
R.K.F. in this case is almost one). 

(2). The difference in gain between the M.K.F. and the 
R.K.F. as a function of -10 log (a) (Fig. 4.2-2). 

The comparison of the gains at TTG is given here. 

As was expected, the gain of the M.K.F. is lower than 
the R.K.F. and reaches the regular gain as a ^0. 

This is an example where the M.K.F. approaches the 
R.K.F. in the limit. 

b. Influence of the Mean of X(0) to the 
Performance of the Filter'" 

The height and the velocity of the chaff are 
taken in this case as parameters. Both are distinguishable 
parameters between the target and the chaff. Thus we 
study in this simulation how they influence the performances 
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Fig* 4.2-2 The difference in gain between M.K.F. and 
R.K.F. as a function of • • 0 Cou, [ a.) ■ 



i 

! 

i 



of the filter. The most distinguishable parameter is the 
initial height. The following parameters and matrices 
are taken for this test: 



1 0 0l 



$ = 



0 

0 



1 1 
0 1 



1 



c = 



0 



0 



1 



0 



0 



1 



R = 



0 



0 



1 



M 



'50 0 0 " 

0 5000 0 

0 0 0.5 



Q = 



1.0 

0 

_ 0 



0 

9-10 

0 



0 

0 

0.25 
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I 



X(0) 



10 “ 

a 

-1 



where "a" is a parameter to be changed. 

Two results are given here: 

(1) The probability of missing a's a function of 
"a" (X2)(0)) (see Fig. 4.2-3). We see 
from the results that the probability of 

a miss is larger than that of the R.K.F. 
without switching, as can be expected, and 
it reaches the performance of R.K.F. in the 
limit (when at the instant TTG the two 
heights are far apart, a > 100 or a < 10) . 
The maximum probability of missing is 
reached at about a ~ 50 (target and chaff 
as approximately at the same height at TTG) . 

(2) The initial gain as a function of "a” in 
Fig. 4.2-4). Here we see that in the limit 
the gains are the same as for the R.K.F., 
and reduce to 0.5 when the two heights are 
indistinguishable as expected. 

The following inputs for testing the influence 



of the velocity are now chosen: 



Probabi 1 i ty of 




100 



Fig- 4.2-3. Probability of miss as a function of the mean 



o 
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Fig. 4.2-4. The initial gain as a function of the 
mean of X- (o) . 



0 



$ = 



1 

0 



0 

1 



0 0 



1 

1 



c = 



1 

0 



0 0 
1 0 



R = 



1 

0 



0 

1 



Q 



1.0 0 0 
0 0 
_ 0 0 



M 



25 0 0 

0 25 0 

_ 0 0 c^_ 



X(0) 



10 

10 

a 



where : 



"a" is a parameter to be changed (the initial 
velocity of the chaff) 
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b 



0.1 a + 0.01 



c = 2b 

Note that in this test the mean of the initial height of 
the target and the chaff are the same. So the only dis- 
tinguishable part is the velocity. The target has velocity 
zero, so the most difficult decision is at velocity zero. 

In order to check the case when the target and the chaff 
have the same dynamics, we choose b and c so that for a = 0, 
the chaff remains at almost zero velocity. For this case 
(chaff velocity about zero) we must expect about 0.5 proba- 
bility of missing. As the difference in the two velocities 
increases, we can expect that the probability of missing will 
approach the probability of missing for the R.K.F. On the 
other hand, for the R.K.F. where there is switching in the input, 
the probability of missing is at the lowest value for 
velocity zero (because both of the targets have the same 
height) , and approaches one as the velocity varies from 
zero. The results f of the simulation are given in Fig. 4.2-5. 
The results agree with the discussion above. 

c. Influence of the Standard Deviation of the 

Initial State to the Performance of the Filter 

The standard deviation of the initial state of 
the chaff changes the probability of the decision at the 
beginning of the process. 

Thus it is interesting to see how it influences 
the performance of the filter. The following matrices 
and parameters are taken for this simulation: 



I 



miss 



O CO V£> ^ CSJ 




Probability of miss as* a function of the mean of X^(0). 



$ = 



1 0 
0 1 
0 0 



0 

1 

1 



c = 



1 

0 



0 

1 



0 

0 



R = 



1 

0 



0 



1 



Q 



1.0 



0 



L 0 



0 0 ~ 
9 * 10 "^ 0 

0 0.25 



X(0) 



10 

100 

- 1 _ 



M = 



50 

0 

_ 0 



0 

0 



0 0.5 



where "a" is the parameter to be changed. 

The difference in probability of missing, 
a function of "a", between the M.K.F. and the R.K.F. 



as 

(without 



switching) is given in Fig. 4.2-6. We see that in the 
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difference in the initial gain between 



limit (a -*■ 0) , the two performances are quite close, and 
when the standard deviation increases, the performances are 
far apart. I.e., the probability of missing for RKF is 
better by 0.1 over the MKF for a = 100). The difference in 
gain between the two filters is given in Fig. 4.2-7. in 
the case 0 ^ 2 ( 0 ) is most influenced by the parameter and 
we show this difference in Fig. 4.2-7. As expected, the 
difference between the gains increases as the standard 
deviation increases. 

d. The Gain as a Function of Time 

To illustrate the behavior of the gain as a 
function of time, a typical example is given here. We know 
that the gains, contrary to the R.K.F. are data dependent; 
so for different data we get different gains. 

The inputs to the system are : 



$ = 



1 

0 

0 



0 

1 

0 



0 

1 

1 



c = 



1 

0 



0 

1 



0 

0 



R = 



1 

0 



0 

1 
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1.0 



0 



0 



Q 



0 
L 0 



9*10 

0 



-4 



0 

0.25 



X(0) 



10 

100 

_-l_ 



M = 



50 0 

0 10 

. 0 0 



0 ■ 
0 

0.5 



An explanation of Fig. 4.2-8 follows. At K = 1, a is about 
0.5, thus the gain is reduced to about .5 of the R.K.F. 
gain. At K = 4 , a becomes almost one and the gain first 
increases and becomes larger than the R.K.F. gain, and then 
it decreases below the R.K.F. gain because a is about 1 
(a = 0.998). It then reaches a steady state gain. At 
K = 35, a is reduced to a low number (a = 0.57) , because of 
the large noise in the input. Hence the gain reduces to 
almost half of the steady state gain. At the next instant 
(K = 36) cc becomes again about 1 and the gain is above the 
R.K.F. gain. After that ct remains almost 1 and the gain 
decreases again and reaches the steady state value. 

For the M.K.F., the steady state value is 
reached only for the intervals where cC is almost one for 
several K's. It can deviate from this value for short 
transients because of spurious noise. 
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3 . Conclusion of the M.K.F. Simulations 

Based upon the simulation results we can conclude 

that: 

a. The performances and the behavior of the M.K.F. 
is as expected. 

b. The M.K.F. is applicable to the chaff-target 
problem and handles the ambiguity problem in 
the monopulse solution (see Chapter IV. A) . 
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V. SOLUTION FOR UNRESOLVED TARGETS FOR C.S. OR L.S. 
- THE PROBABILITY FILTER 



A. THEORETICAL SOLUTION 

In case of L.S. or C.S. it can be assumed that a phase 
detector is not available. Thus we get only amplitude 
from the receiver. The approximated p.d.f. of the return 
signal exists in the literature [ 8 ] , [ 9 ] . Our purpose is 

to approximate the amplitude of the target or the chaff as 
received at each antenna so that the position of the target 
can be calculated. In case of C.S. we can separate the 
signals received during one conical scan period as equiva- 
lent to n separate antennas where n can be as large as the 
number of pulses in one rotation. 

The information of the amplitude of the target at each 
antenna is sufficient to determine the position of the 
target with respect to the antenna (see Chapter II. B). 

The probability filter proposed here seems to be the 
most practical solution to our problem because: 

(1) It gives a good estimation for parameters 
(signal amplitudes in this application) . 

(2) It is easy to implement. 

(3) The probability filter does not depend upon an 
exact knowledge of the p.d.f. of the signals, 
and an approximation is sufficient- This filter 
can be implemented with a premeasured reference 
signal which is used with an approximate p.d.f. 
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to obtain the solution. This is an advantage over 
a likelihood filter which cannot be implemented here 
because the p.d.f. of the signals is not know exactly 
The theory for the probability filter is developed in 
Appendix D. The application of this general approach follows 
with simulated examples to demonstrate its effectiveness. 
Because we restrict ourselves to surface targets we can 
assume a stationary process for the amplitudes measured at 
the receivers. 

Thus we can use (D-12b) as an approximation to the time 
average (i.e. with reasonable "T" compared to the missile's 
flight time) . 

B. PROBABILITY FILTER SIMULATION 
1. Introduction 

In this section we give some results using the proba- 
bility filter concept. The difference between the likelihood 
filter and the probability filter is that for the likelihood 
filter we need more information about the process, i.e., we 
need to know the joint p.d.f. ( f ^ (x^^ , X 2 r • • • / x^) , while in 
the proposed probability filter we have to know only the 

marginal p.d.f. (f^ (x) ) or its approximations. The relation 

^i 

ship between the samples involved the expected value of a 
function of the samples, and for this the sample average 
is used (assuming an ergodic process), i.e.: 

E{f(Z)} = f(Z) = I Z f(Z(i)) (5.2-1) 

where the symbol means time average. 
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In order to estimate a parameter vector 0 we develop 
the theory of the probability filter given in Appendix D. 

The main idea of the probability filter is to minimize the 
difference/ in some sense, between the conditional p.d.f. 
f[Z|0] (where 0 is the stimate of 0), and the measured p.d.f. 
of the process, Z, which contains the information of 0 , 
f(Z) = f(zl0). We choose the norm to be (D-11) : 

00 

I = j j[f^{z\B) - dz min. 



This leads to a necessary condition (D-12) : 



“3f„(z|§) ^ ”9f_(Z|0) 

° = - / - ^-3e- f.(z)dz 

90 ~ 

—00 —00 



Our task is to solve this equation for 0. 

Summary of notations : 

0 = estimation of 0 

0 = true value 

9 j _ jth iteration in the process of finding 0. 

0 (K) = 0 at the end of N data block. 

/N /V 

From D-12 we want to determine 0 = 0 so that 31/30^ = 0 
for all i. In order to accomplish this the following procedure 
is used; (D-11) and (D-12) can be rewritten as 



1=1 / [f2(Z|Y) - f2(Z)]^ dZ 



(5.2-la) 



, oo3f (ZlY) -3f (ZlY)^ ^ 

Iy = / -A.„ ..:^..f^(z|Y)dZ - / ,^~- f,(Z)dZ ^ f^(Y,0) 



3Y 



3Y 
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The experimental value of 3I/3Y is given by 



[— ] = f 

E J 



”3f^(ZlY) N 3f„(Z(i) lY) . 

— a -- - 



i=l 






(5.2-lb) 

Note that fg(Y,0) does not depend explicitly on 9, although 
9 of course influences the Z(i). 

/\ 

We are seeking 9 such that fg(9,9) = 0. Usually we 
solve this iteratively; 0^ is guessed, and then is 

/N ^ 

obtained from 0.,, = 0- - £’f„(0.,0) for j > 0, where e is 

/N /% 

an arbitrary number. If the sequence 0j converges to 0, 
then clearly fg(0,0) = 0. If z is very small, the convergence 
is slow. On the other hand, if £ is large we might get 
oscillations. Thus we must use for this method an optimal 
s, which might be a function of f_. A different way to 

^ II, 

solve for 0 is to use Eqs . (5.2-la) and (5.2-lb) that can 

be iterated so that 






(5.2-lc) 



where 0. is the j approximation for 9. Eq. (5.2-lc) is 

/N 

solved successively for ri-1 steps, until 

^ ~ '' A ~ 

f„(0 ,0) = 0, and we conclude that 9 = 0„ = 9. 

After we iterate (5.2-lc) for the first N samples of Z(i) 

we get 0(1). Then we take a new sequence of N samples and 

/N 

(5.2-lc) is iterated again to produce 9(2). Thus after the 
probability filter we get a sequence of the estimated parameters: 



{0 (1) , 9(2), 



0 (K) , . . . } . 



• • • / 



If actually 0 (K) changes between blocks of N samples, we may 

/N 

want to smooth the estimates 0 (K) with a Kalman filter. In 
the discussion which follows the following terms are used: 
"Open loop" means the estimation of 9 (K) with a 
single block of N samples Z(i). 

"Closed loop" means the estimation of 0 (K) with 
the "open loop" followed by a Kalman filter. 

As we know (see Appendix F) , the outputs from 
the Kalman filter are 

0(k|k) - the update estimation of 0 (K) 

0(k|k- 1) - the predicted estimation of 0 (K) 

The concept of the "closed loop" is discussed and 
simulated for the estimation of one parameter (parts 2 and 3) . 
2 . Open Loop Performance - One Parameter Estimation 
In this section we give the open loop output for 
estimating one parameter. We choose the exponential 
distribution. The p.d.f. is given by; 

f_(Z|0) = 0 e"®^ Z > 0 (5.2-2) 

0 is the parameter to be estimated. The equation for the 
probability filter is as follows: 



9f2(Z| Y) 
3Y 



-YZ 

e 



YZ e 



-YZ 



(5.4-2a) 



c°3f (z|y) 

J ■ - V ' 



/ Y e”^^^dZ - / Y^Ze ^^^dZ = 0.25 

0 0 

(5.2-2b) 



Thus we can rewrite equation (5.2-lb) for our case by 
substituting (5.2-2a) and (5.2-2b) into (5.2-lb). 



i 

I 



(5.2-3) 



N 



f^(Y,9) = [||] = 0.25 I [l-YZ(i)]exp[-YZ(i)] 



N i=l 



The first term in (5.2-la) is given by: 






8Y 



dZ 



0 



0 



= [9/(9+Y)] 



(5.2-4) 



Note that (5.2-2b) is a special case of (5.2-4). From 
equations (5.2-2b) and (5.2-4), (5.2-la) can be written 

for our case as 



f^(Y,6) ^ [|i] 



= 0.25 - [9/(9+Y) ] 



(5.2-5) 



Now we iterate our solution according to (5.2-lc) using 
Eqs. (5.2-3) to (5.2-5). This leads to the solution: 



9j^^ = 9j [Uj/d-a^) ] (5.2-6) 



where : 



a 



- ~\l- 1 



[l-9jZ(i) ]exp[-9^Z(i) ] 
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Note: If the argument under the root turns negative then 

one may take a larger N. 

To show that this concept works, we now test the "open 
loop" criterion. For this simulation 9 = 1 is taken. We 
perform three simulations; 

- Testing 31/ 3Y 

The purpose of this simulation is to show that 

the experimental 3I/3Y contains the information 

required for the iteration. 

Estimating 9 in an open loop. In this simulation 

9j is iterated and we try to find; 

(i) The convergence of the iteration by 
(5.2-6) . 

(ii) The number of iterations needed for 
estimating 9. 

(iii) The variance and the bias error of 

/N 

9. 

Estimating of 9 in a closed loop. 

The open loop estimate of 9 acts like a measuring 
device for 9. The R.K.F. is used here as a filter for 
this data and we test the performances of the overall 
system. The design of the R.K.F. is influenced by 
the variance of the error and the bias in the output 
from the open loop probability filter, i.e., we assume 
that after the probability filter we have a laeasured 
signal with noise. 
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9 (K) + V(K) 



(5.2-7a) 



9 (K) = 



where 

9 (K) is the true value, and 

V(K) is the measurement noise at instant K. 

In this case the statistics of V(K) are a function of 
9 (K) , and this function can be evaluted approximately by 
assuming that 9 (K) = 9 (K) . We are interested in the 
variance of V(K) which is denoted by R(K) . Thus, 

A 

R(K) = VAR[V(K)] = R[9(K)] := R[9(K)] (5.2-7b) 

The process 9 (K) is assumed to be described by a linear 
difference equation similar to (F-2) and the statistics 
similar to (F-3,4,5). I.e., 



9 (K+1) = ^9 (K) + W(K) 

/V 

9 (K) = 9 (K) + V(K) 



(5.2-7C) 



where, in (F-2), 9 (K) replaces X(K), and 9 (K) replaces Y(K). 
The R.K.F. can be implemnted only when 

R(K) = R[9(K)1 = Rl.9^(K)] , 

M, Q, and 9(0) are given. 
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a. Testing the Experimental 3 i/9y 

The purpose of this simulation is to show that 
3I/9Y contains sufficient information to close the loop by 
iteration. 

In this case we take the following parameters : 

(1) 9 = 1.0 

(2) Sample size - 1000 

(3) Number of ensemble members = 1. The 
results are given in two figures; 

Fig. 5.2-1 and Fig. 5.2-2 which give the 
results of 9 i/3y for 9^/9 = 0.5 to 5.0. 

It is significant to note that; 

/N 

(i) At 9^/6 = 1 the error between the 
experimental gradient and the analytic 
gradient is almost zero. 

(ii) The values of the results are very close 

/N 

to the analytic values for all 0 . 

o 

3 1 

(iii) The curve for approaches the 

asymptotic line 0.25 which is the 
correct analytic value. 

The results show that the concept of the 
probability filter works. 

b. Estimation of 0^ in an Open Loop 

From the first simulation we have seen that 
9 i/3y contains information about 0. The purpose of this 
experiment is to study the bias error and the variance of 
estimation of 0 for two cases ; 
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rigure 5.2-1. The open loop output. 
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Figure 5.2-2. The error li> the open loop output. 



when 9^ = 0 we will see the influence of 

the number of samples on estimating 0, by 

/\ 

one iteration (9^). 

/N 

when 9^ is a parameter, we will see that 
the open loop estimator always converges 
toward the true value of 9 . 

In both of these cases the number of ensemble members is 
500. The results are given in Figures 5.2-3 and 5.2-4 
and in Tables 5.2-1 and 5.2-2. 

Refer to the first simulation, i.e. one changes 
the number of samples, N, of (5.2-3). According to the 
figures it seems that the estimator behaves similar to the 
likelihood estimator, i.e.: 

- The estimator seems asymptotically unbiased 
(see Fig. 5 . 2-3) . 

The variance of the error reduces inversely 
with N (see Fig. 5.2-4). 

Refer to the 2^^*^ simulation, i.e., we change 

/V 

9^ as a parameter, and look one step forward to see the 

/N 

error of 9^^. We define a new variable: 

a(j) = 9j/9, j = 0,1,2, ... (5.2-8) 

The number of samples, N, for this simulation = 1000. 

Refer to Table 5.2-1. From the table we see 
that the open loop estimator gives a very good estimate of 
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Fig. 5.2-3. The bias error as a function of the number of samples. 
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TABLE 5.2-1 



THE ERROR AND VARIANCE OF ESTIMATING 

/N /S 

0, WHEN 0 IS A PARAMETER 

1 O 



/V /V /V 



a(0) 


eU.o 


- a(l)} 


VAR{a(l) } 


0.1 


-1.2- 


1 

o 

f— 1 


9.4 


• lo"^ 


0.2 


-1.1 


•lo"^ 


1.1 


• 10“^ 


0.5 


-1.9 


1 

O 
• — 1 

• 


1.4 


• 10"^ 


1.0 


-2.0 


• 10 ^ 


2.1 


• lo"^ 


1.5 


1.7 


ro 

1 

o 
1 — 1 

• 


3.1 


• lo"^ 


2.0 


8.7 


1 

o 

f— 1 

• 


4.2 


ro 

1 

o 
1 — 1 

• 


5.0 


1.8 


♦ lo"^ 


2,2 


• lo"^ 
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TABLE 5.2-2 



0 

0.1 

0.2 

0.5 

1.0 

2.0 

5.0 

10.0 



THE INFLUENCE OF THE ABSOLUTE VALUE 
TO THE RELATIVE ERROR AND VARIANCE 



Absolute Value 



Relative Value 



E{0-0,} Var{0,} E{l-a(l)} Var(0,/0 ) 



-2 • lO""^ 2.06*10"^ -2 • lO"^ 



2.06-10 



-3 



-4.5-10"'^ 7.5 -lO"^ 2.3 • lo"^ 1.88-10 ^ 



2.9 -lO”"^ 5.2 -10“^^ 5.8 * lO”"^ 2.06-10~^ 



-2 • lo"^ 2.06*10“^ -2 ♦ lO”^ 2.06*10"^ 



9 • lO""^ 7.9 -lO"^ 4.5 *10~^ 1.96-10"^ 



5.6-10"^ 5.1 -lO"^ 1.1*10"^ 2.05-10"-^ 



5.9-10 



-3 



0.206 



S.O-IO""^ 2.06-10"^ 
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the parameter, even if we are far away from the true 
parameter value initially. Thus it seems that a few itera- 
tions are good enough to estimate the parameter since the 
single iteration used here gives such accurate results. We 
can use two criteria for determining the number of iterations : 
a fixed number of iterations or a flexible number of itera- 
tions. For example, a stopping criterion might be 



/N 

jl.O - a*(K) I £ e 



(5.2-9) 



where 

a*(K) = (5.2-10) 

®K-1 

« 

In this example, it seems that a fixed number of iterations 
are enough. A simulation is now performed to prove this. 

/V 

We start with 9 =9 and calculate 9, for a 

o i 

range of values of 9. The results are given in Table 5.2-2 

2 

and indicated that the relative variance of 9^/9 is 
independnet of the value of 9 . The variance is about 

-3 

2*10 for all the tested values. And the relative error 

-3 

IS a very low mamber and it is of the order of 10 . Thus 

we can conclude that our results are generally independent 
of the values of 9 . 

To prove that a fixed number of iterations are 
enough (two iterations give a very good estimate) , we make 
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several runs with different initial conditions. The nuinber of 
iterations is 10 andN=1000 . The results are given in Table 5.2-3. 

From the table we see that there is no significant improve- 
ment after two iterations, even if our first value is far 

/\ /N 

away from the real one (a(0) = 100 or a(0) = 0.01). Thus 
for the estimator we use two iterations. From the table 
we see that we get a bias error on the order of 0.2%. 

The relative variance is given by: 

VAR[a(2)] = 0.002 

Thus we can say that 



VAREe^] = 0.002 0^ 



and can be approximated by; 



/N /N 

VAR[02] = 0.002 02 (5.2-11) 

This equation will be used for the noise in the closed 
loop simulation which follows. 

3 . Closed Loop Performances - One Parameter 

After we have examined the open loop performance 
of the filter, we can close the loop with a Kalman filter 
and estimate the parameter as a function of time. The 
conceptual block diagram is given in Fig. 5.2-5. The con- 
cept of the filter is very simple. The open loop estimate 
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Test of the Improvement by the number of iterations to the estimation 
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Figure 5.2-5. The close loop estimation of 0 



of 9, acts as the input to a R.K.F. The difference equation 
of 9 is given by (5.2-7c). The random measurement noise is 
data dependent according to the simulation in the last para- 
graph (5.2-11). So the variance of the estimation noise, 

R(K) , is computed on line. Hence the Kalman gain must be 
computed on line. 

Explanation to the figure: We use the open loop 

estimator for 9^ as a measurement of plus noise. The 
initial guess for 9, Q 2 ' the predicted state from R.K.F. 

/V 

(i.e. 9 q = 9(k|k- 1)). We compute on line the variance 
according to the assumed R(K) using equation (5.2-11); 

R(K) = 0.002 9^^ (5.2-12) 

The R.K.F. receives the input 9^, the computed R(K) , and we 
then estimate the state 9(K|K) and the predicted state 9(KlK-l). 
In summary; N = 1000 samples of the data are taken, two 
iterations are used to estimate 9(K), so that the Kalman filter 
updates 9(K) in probability filter in successive blocks of 
1000 samples. The random generators which provide W(K) and 
9(0) are Gaussian (see 5.2-7c). Since 9 (K) is positive the 
above procedure is made in order to avoid this problem: 

r b if b > 0 

6(0) = \ _ 

I 9 (0)/5 if b < 0 



where b is the output from the generator which produces 9(0). 
And ; 



9 (K+1) 



(|)9 (K) + W(K) if it became positive 
9(K)/5 elsewhere. 
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The first simulation that we perform is a case where we 
have a process which has the following parameters; 

$ = 1.0 

Q = 0.1 

■ 0 ( 0 ) = 1.0 

VAR{0(O)} = 1.0 

For the Monte Carlo simulation 250 members of the 
ensemble are taken. The performances come out as expected 
with: 

- Variance of the error over the ensemble is of 
the order of 4 • 10 for all K. 

- Bias error over the ensemble is of the order of 
2 • 10~^, for all K. 

The second simulation is identical to the first one but 
with different parameters. In this simulation we want to 
examine two things : 

- The estimation error bias and variance as a 
function of K. 

- The theoretical variance of the random measure noise 
which is computed according to the experiment formula 
(5.2-12) by comparing the calculated variance in the 
estimation using (5.2-12) and (F-10) with the measured 
variance of the estimation of an ensemble of 250 
experiments . 

The input parameters are; 

$ = 0.9 

Q = 0.5 
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9(0) 



5.0 



VAR (0(0)} = 2.0 

Ensemble members = 250 

The results (see Fig. 5.2-6 and Fig. 5.2-7), show 

that; 

- The variance of the estimation error is on the 
order of 2 • 10 

-2 

- The bias error is on the order of 10 

- The theoretical variance of the estimated error 
is very close to the simulated variance. Thus 
our experimental formula (5.2-12) is accurate. 

- We know that our error is biased but it 

is one order less than the standard deviation of 

the error. Hence the bias error is negligible. 

4 . Estimation of Several Parameters by the 
Probability Filter 

a. Introduction 

In the last two sections we estimated one 
parameter in an open loop and a closed loop simulation. In 

this section the open loop probability filter performance is 
tested for more than one parameter. For multiple parameters 
it appears that we often cannot find a close form formula 
for iteration procedure for the parameter values, some have 
to solve for them numerically. Performances of the probability 
filter is compared to the likelihood filter. 
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theoretical variances of the estimation error 
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b. Mathematical Equations 

As an example we take a Gaussian p.d.f. with 
two unknown parameters 

P - the mean 

2 A 1 

“ ~ 2 

a 



where: 



o - is the variance of the process. 
The p.d.f. of this process is given by 



f^ (Z/a,\i) 



exp{- 



(5.2-13) 



Let's define the parameters*; 




(5.2-14) 



Define a new variable related to the real 

/N /N 

parameter (a and y) and the estimated parameters ( 0 j^, 02 ): 
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ft2 A Q 2 2 

3 = 02 + ct 



A 1 ^ 2^ 2 

Y = ^[9-^9, + a^y] 

6 



^2 A ra 2 ft 2 ^ 2 2, 

c = [02 9^ + a y ] 



(5.2-15) 



K A 2 2 ^2 

b = c - Y * B 



9. 



2 2 



d = 1 - 9„ Y + 2 y9,9., - 0, 0 



2 : 2 



B 



12 1 "2 



9 /V .. 

= 1 - -f- - 0 /(y-0t)^ 

3 



The equation of the probability filters follows 
according to (5.2-la,b,c) as: 



3f, 

w. 



— (Z-Y, )exp{- iY_^(Z-Y^)^} 
/2^ ^ 2 2 1 



3f, 

W. 



— [1 - Y 2 ^(Z-Y^) ^]exp{- -|y 2 ^ (Z-Y^) 
/2tt 



Thus from (5.2-lb) we get for this case: 
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[^] 

E 



[3I_] 



^ N 3f2^(Z(i) |Y^,Y2) 

~ N" . ^ Y 

i=l 1 



, N 8f„ (Z(i) |Y^ ,Y.) 

1 - i y _ Z '12 

4/? 



3Y, 



(5.2-16) 



The theoretical 3I/3Y^ comes out to be. 



according to (5.2-la): 



[1^1 






[^] 

n, 

2 T 



-i^fexpt- |bl-[Y - V^l 

/ 2tt 



a 



_ 4 /tt /2Tr • S 



exp [ - ] ’ d 



(5.2-17) 



We 



cannot find an iterative formula relating ^j + 1 



to ttj and Uj directly. However, we still need to be able 
to solve the relationship 

/N. /N. 

x'N 



ft(0.9) = 



~hj ~ ~ 



(5.2-18) 

- ^ ^ 

0=0 

Assume that we are close to the real parameters and then: 



/V 




U 
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1 



©2 = a 

Y = M 

b = 0 

a 

3 

d 

We estimate the parameters in the experimental 
equation as suggested in Appendix D by using a reference 
signal (see Fig. D-5) . 

Our system measures [8l/8Y-]„. The analytical functions 
[8I/8Y^]^ are given by equation (5.2-17). Using the form of 
(5.2-lc), (5.2-18) leads to the approximate equations 







'"j+1 






- [ 



31 



2 /tt 
a-i 



a 



2 

j+i 



a 2(1 _ [ 1 ^] . 4 /?) 



(5.2-20) 



We see that in the limit, i.e., when [3I/3Y-]„ = 0, we 

j- hi 

get the estimation of the parameters . 

It can be seen from (5.2-17) that the sign of 

the error in t is independent of , and it is dependent 

/\ 

only on the error (y - u) itself. Consequently, the 
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convergence is better if is fixed and is changed until 

/\ ^ 

[3I/3y. ]„ ->■ 0. Then y. and a. are changed until [3I/9Y,] 

J J D IE 

and [3I/3Y2]g approach zero. In this case it appears that 
we can not use a fixed number of iterations, we have to 
iterate until [3I/3Y. ]„ ^ 0. 

X XI 

c. Open Loop Performances 

In the last paragraph it has been shown that it 
is sufficient to treat the open loop probability filter as a 
measurement device and add a Kalman filter to it with measure- 
ment noise R(K) dependent on the data. We perform two experi- 
ments, in which we use the probability filter with a reference 
signal (see fig. (5.2-9,)). 

The first experiment is to show that the performance 

of the filter is independent of the mean, U* The second 

-2 

experiment tests the influence of the variance, 6^ on the 
perfomance of the filter. 

In each of these cases, our ensemble number is 50, 
and N = 1000. For reference we include the performance of 
a likelihood filter for competition. The results are given 
in Tables (5.2-4) and (5.2-5). 

The first experiment (see Table 5.2-4) shows clearly 
that the error in the variance of 9^^ and ©2 is essentially 
independent of 9^^. We hold ©2 fixed (in our case ©2 = 1.0) 
and changed ©^ from -50 to +50. The results for each value 
of ©^ are very close. 

The parameter ©^^ is chosen to be zero and ©2 is 
changed from 0.01 to 100. From Table 5.2-5 we see that the 
variances can be expressed empirically. 
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/N 

Fig. 5.2-?. Block Diagram of estimating 9. 
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Definitions for the folloving tables; 
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Table 5.4-5. Tbe influence of 02 on estimating 0, and ©2 should be 0 in theory 



- p 

VAR ( 0 ) 



= 1.4 • lO"^/02^ 



(5.2-21) 



VAR(02^) = 1.3 • lO"^ 02^ 



If we compare the results of the probability filter to the 
likelihood filter (Table 5. 2-4,5) we see that the performances 
are close, but the likelihood filter gives better results. 

The experimental results, which are given by 
formula (5.2-21), must be used when closing the loop with 
a Kalman filter as discussed in the last section. 

We will now examine the influence of the momber 
of samples, N, on the estimation for open loop probability 
filter. The number of samples has great influence on the 
estimation process because as the number of samples increase, 
the accuracy of the estimation increases. But, we have to 
look at the whole system, i.e.: 

The dynamics of the process generating Z(t) 
when the dynamics of the process are rapid, 
we would like to take a minimun number of 
samples . 

A good estimation after N samples does not 
mean that we cannot get a better estimation 
with fewer than N samples (for example 
N/10) . If we can filter the output of the 
probability filter by a Kalman-type filter 
we may get a better estimation. For example 
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N = 5 might be sometimes better than N = 1000 
if we place a Kalman-type filter after the 
probability filter. 

Thus, the information about the influence of the number of 
samples is very important when we are going to design a 
whole system. The results of the test of the influence of 
the number of samples, are given in Table 5.2-6 and Figures 
5.2-9 and 5.2-10. From the results we see that the variances 
of the estimation errors decrease almost inversely (except 
for N = 2) v/ith the number of samples. From the results it 
seems that the variances of the probability and likelihood 
filters decrease almost inversely with the number of samples. 
The likelihood filter (except at very small N) gives results 
for most of the cases better than the probability filter. 

d. Estimation of Parameter in Case of p.d.f. 

Not Known Exactly 

As was discussed in Appendix D, by having a refer- 
ence signal we can estimate parameters even if we have only 
an approximation for the p.d.f. (See equation D-12b.) 

In this section we would like to show how this 
concept works. We take at first a Gaussian process including 
a limiter, which is very common in any practical measurement. 
We want to estimate the parameters with a minimum influence 
of the accuracy of the assumed p.d.f. We test this case with 
the probability filter with reference signals, and we compare 
it to the likelihood filter. The block diagram of the pro- 
cess is given in Fig. 5.2-11. 
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Table 5.2-6 The influence of the number of samples (N)' to 
the performances of the filters 



I 

( 




Fig. 5.2-9 . 



The variance of 6, as a function of the 
sample's number. ^ 
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Fig. 5.2-10. The variance of 02 a function of the 
sample ' s number . 
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The block diagram of the process with limiter 



For the probability filter the reference signal 
is taken to be a limited Gaussian signal and the assumed p.d.f. 
is Gaussian. For likelihood filter the assumed p.d.f. is 
taken to be Gaussian also. 

We test first the case of a symmetrical limiter, 
i.e., |Z . I = Z . The mean of the output is not influ- 
enced by the limiter. The estimation of 02 is influenced by 
the limiter. The bias error of 8 ^ strongly influenced 
by the limiter. The following parameters are taken for this 
test . 



# of samples (N) = 100 ; 

# of ensemble members = 100 ; 

From the first simulation (see Fig. 5.2-12) we can conclude 
that : 

(i) The probability filter is much better than 

the likelihood filter when Z <3; for 

max 

Z >3 the likelihood filter becomes a 
max — 

better filter than the probability filter, 
(ii) The bias error of the likelihood filter 

depends very much on the knowledge of the 
p.d.f., while for the probability filter 
error it is almost constant. 

The second simulation (with a non symmetrical limiter is 
performed), i.e., I Z . I Z . We choose Z . = -1, and 

we change The unsymmetrical limiter influences the 
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Os. 



10 . 
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A0^ = The bias error of estimating 
fP-Probabi 1 i ty filter (• ) 

Ll -b,i kl i hood filter (O) 

Z =-Z . 
max min 



O 



« 







10 -I , 

0.01 0,1 



z 

ma 



10 . 



Fig. 5.2-12. The bias error of ©2 in 

case of symmetrical limiter. 
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TABLE 5.2-7 



Z 

max 

0 

0.01 

0.5 

1 

5 

10 



THE INFLUENCE OF NON- SYMMETRICAL LIMIT 
ON ESTIMATING THE PARAMETERS 



A9^P 



Ae/ 



A9, 



A9, 



5.5 • lO"^ -1.3 • 10"^ 0.32 



-1.5 



5.6 • lO"^ -1.4 • 10"^ 0.31 



-1.5 



1.1 • lO"^ -1.3 • lO"^ 0.11 -0.74 



3.8 • lO"^ 2.4 • lO"^ -2 • lO"^ -0.41 



1.4 • 10“^ -4.1 • lO"^ -7.2 • 10“^ -0.18 



-6.7 • 10 ^ -7.6 • lO"^ -9.0 • 10 ^ -0.16 



Z . 
min 



-1 
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variances and the means, so we have to check the bias error 
of estimating 0^ and 02. The results are given in Table 
7.2-7. The results show that the bias errors of the proba- 
bility filter are influenced weakly by the p.d.f. while for 
the likelihood filter there is a great dependence upon the 
knowledge of the p.d.f. 

It should be mentioned that the likelihood filter used 
for comparison here is a "straw man" in the sense that the 
probability filter is given the correct class of p.d.f 's 
(indirectly through the reference signal) , whereas the like- 
lihood filter is not. In this example, one could design a 
likelihood filter using the correct class of p.d.f. 's that 
would probably outperform the probability filter. However, 
our point is that the probability filter is robust (in this 
example - further work is in order) with respect to errors 
in establishing the class of p.d.f. 's used for gradient 
computations when a reference signal is available. In prac- 
tice, the choice between filter types would depend on whether 
the correct class of p.d.f. 's can be well approximated 
(likelihood filter) or not (probability filter) by simple 
analytic expressions. 

5 . Conclusion of the Probability Filter Simulations 

After the simulations we conclude the following: 

a. It has been proven that the probability filter 
can be implemented quite easily. 

b. The performance of the probability filter is 
(in our example) not sensitive to errors in 
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estimating the form of the p.d.f. f„(z|0) when 

Lk ~ 

a reference signal is available. 

c. The performances of the filter when the p.d.f. 
is known exactly, is of the same order of magni- 
tude compared to the likelihood filter, but 
gives worse performance most of the time. 

d. In the target/chaff problem we do not know the 
exact p.d.f. because there are several types of 
chaff and targets. It is quite easy to provide 
a reference signal for the probability filter 
which will fit to the exact p.d.f. Thus the 
probability filter may be found to be the 
practical solution. Further investigation is 
needed . 
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VI. SUMMARY AND CONCLUSION 



In this research we tried to solve a realistic and 
practical problem of separation of time signals when they 
have almost the same power spectrum. As an example, we 
took the chaff-target problem in tracking radar where Doppler 
separation is impossible. We started with a discussion of 
the background of the problem and developed two possible 
solutions . 

a. The M.K.F. with antenna modification. This is one 
solution to the problem in the case of maximum 
information available (for example in the case of 
monopulse radar with phase detector) . 

b. The probability filter. This is one possibility 

of the solution to the problem in case that minimum 
information is given (in the case of concical scan 
radar or lobe switching) . 

To demonstrate our solution to the problem we have proven 
by simulations that the two main filters, the M.K.F. and the 
probability filters work as expected. 

The following items remain to be investigated as an 
extension of the work. 

a. The complete solution in case of C.S. and L.S. 

b. The optimum number of samples (N) for the 
probability filter as a function of the dynamics 
of the signal process. 
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c. The probability filter extended to nonergodic 



processes . 
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APPENDIX A 



MODIFIED KALMAN FILTER FOR THE CASE THAT THE MEASUREMENT 
COMES FROM TWO DIFFERENT SOURCES 

1. Introduction 

In many cases we have to filter data which is not 
the regular data assumed by the Kalman filter. The Kalman 
filter is well established and sometimes it is convenient 
to modify the problem such that it will fit the Kalman 
filter, so we will have a modified Kalman filter. In our 
case (see Chapter IV) such a problem arises. The problem 
can be described briefly as follows (refer to Fig. A-1) : 
we have a regular process described by (F-2) which gives 
us the output vector Y(K). Y(K) can be split into two 
vectors (Y^^^(K) and Y^^^(K)): 

Y(K) = 

where Y^^^ and y^^^ are vectors of the same dimension. 

The vectors are transfered through a random process, 
which switches between them. After the switching block we 
get our measurement vector, Z (K) . The vector Z (K) can be 
equal to one of the two possibilities: 



Y^^^ (K) 
( 2 ) 

(K) 



Z (K) 




(K) 



(K) 
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)(K) 




j, 
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Figure A-1. The measurement process. 



or : 



Z(K) 



(K) 



(K) 



The problem statement is: given the measurement Z(l), .../ 

Z(n) find the state vector X(n). 

Four papers deal with a similar problem [26,27,28,29]. 
The first two papers cannot be implemented straightforward 
because both of them formulate the problem in a different 
way. The first one [26] deals with multi-target measure- 
ments giving information about the location of each target. 

The system must identify the type of the target and estimate 
its position which is a more general problem than ours, 
because we assiime that we know the nature of the targets. 

The second and third references [27,28] deal with 
tracking in a cluttered environment. This could be adapted 
for our problem (by treating the clutter as chaff) but the 
assumption is that the undesirable returns occur completely 
at random which is not the case in our problem. There is 
no underlying dynamical process from which the returns are 
generated. Hence, no prediction from past data can be made 
on the location, nature or number of these returns at the 
next measurement time. 

The last reference [29 ] is close to our problem but 
uses a different starting point. It does not give the 
optimum filter, which is known to be: 
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oo 



X(K|K) = / Xf(XlZ^^^) dX 

— 00 



(for minimum variance estimates, see Appendix F) . The 

( K ) 

vector Z is all the measured data up to time K. The 
p.d.f.'s f(x|z^ ) are not Gaussian, which makes the evalua- 
tion of X(k|k) very difficult to implement. It is extremely 
difficult to implement on a computer, since the p.d.f. is 
conditioned on all available measurement data. Thus [29 ] 
uses approximate p.d.f.'s. Our goal is to develop a filter 
similar to the Kalman. In order to achieve this we use: 

A selection block which transforms the measured data 
for input to a Kalman like filter. 

- The assumption that all the p.d.f.'s are Gaussian 

i 

(like in [29 ] ) .* 

Refer to Fig. a- 2: we try to reverse the switching 

in the data in the first block, by selecting one from the 
two possibilities of Z(K). Thus after this block if our 
selection was correct we have the correct input to the 
Kalman filter. In a regular Kalman filter (R.K.F.) the 
input components (Z^^^ (K) and Z^^^ (K) ) are not switched. 

Hence in our case if we made an error in the selection the 
output is biased. 

In our problem we have a positive probability to 
select a wrong combination of Z^^^ (K) and Z^^^ (K) . We 
shall modify the R.K.F. so that it will be able to handle 
this error. The next block is a modified Kalman filter (M.K.F.) . 
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This leads to a form like a Kalman filter but with 
modified gain and covariance matrices. In the R.K.F. we 
have two types of estimated states. 

* I 

One IS X(K|K) which is called the update estimated 
state. This is the optimal estimation of X(K) taking 
into account all the observations so far (i.e., at times 

up to and including K) . 

The second estimated vector is X(k|k- 1), which is the 
optimal estimation of the states given all the 
observations up to the last measurement (i.e., the 
observations at times up to and including K-1) . 

This vector is called the predicted estimated states. 

The M.K.F. has similar types of estimated states, 

an update estimated state and a predicted estimated state. 

< 

Intuitively we can say that: 

If we trust our new measurement, we can use the 
update estimated state. 

If we do not trust our new measurement, we have to 
use the predicted estimated state. 

If it is in between, the output must be weighted 
somehow between the two estimated states. 

The weighting depends on the probability of a 
correct solution. 

The estimated states are not optimal, in the sense 
of minimum variance of error, because of our 
assumption that the p.d.f.'s are Gaussian. 
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2. Modified Kalman Filter 



Assumption 

(i) We have two separate, independent sources of 
states with a given dynamic: 

(K+1) = (K) + A^^^U(K) + W^^^ (K) 

^ 

(A-1) 

X^^^ (K+1) = (K) + A^^^U(K) + W^^^ (K) 



(ii) The measurements are given by 



Z (K) = 



L(K) ] 



L(K) = < 



I 0 

0 I 

0 I 

1 0 



^ L. 



^ L, 



X^^^ (K) 



X^^^ (K) 



if H. 



+ V(K) 

(A-la) 



(A- lb) 



if H, 



(iii) The statistics of the random variables involved 
here are (assume all are Gaussian processes) : 



E[V(K)] = 0 4 K 



E[V(K)V^(j)] = R(K) 6 ^. ^ K^j 



(A- 2) 



W^^^ (K) 



W^^^ (K) 



= 0 V K 
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E[W(K) -W (K) ] = 



(K) 1 0 



I (2) 

I ’ (k) 






= 2(K)«Kj ^ ■',3 



( 0 ) 



( 0 ) 



X 

~o 



COVAR 



X^^^ (0) 



X^^^ (0) 



= M 



E[X(0)-V (K)] = 0 Y K 



E [X(0) -W" (K) ] = 0 y K > 0 



Prob (Hj^) 



Prob (H 2 ) = 2 



No correlation exists between H(k) at time K to H(j) at 
time j for all K j . In other words, there is no 
a priori information about the switching between 
and H 2 . 
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{ 



For convenience we define: 



X(K) ^ 



(K) 



(K) 



r,(i) 



$ = 



0 



* ( 2 ) 

0 1 <)) J 



A ^ 



r,(i) 



( 2 ) 



W(K) = 



(K) 



( 2 ) 

(K) 



(A-2a) 



Solution 

We would like to have a sequential estimator 
because we like to estimate in real time with minimum delay. 
The non-sequential estimator (Bayesian approach) will be 
more accurate (i.e., with minimum variance error) but much 
more complicated, needs more storage, and has delayed output. 
For these reasons we will restrict ourselves to sequential 
estimators. Also, we will use a linear estimator. Because 



I 



164 



I 



I 

i 



the Kalman-like filter is linear, and we have to store only 
the last information about our estimation, we will build 
our estimator in two steps - the first one will be a selec- 
tion block and the second one will be a modified Kalman 
filter. Because the Kalman filter is optimal for second 
order statistics and is a linear estimator, the only 
thing that we have to store is the previous values of the 
(i) Covariance matrices and (ii) Predicted state. 

The modified Kalman filter has two main blocks 
(see Fig. A-2) : 

a. The Selection Block. 

In this block we make the decision H = or 
H 2 . Let's define 

y(l) 4 

(A-3) 

^( 2 ) 4 ^i2) ^( 2 ) 



By this definition we can write: 



r 



Z (K) = 



L Y 



( 2 ) 



t-^(2) 



( 1 ) 



if is true 



if H 2 is true 
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Fig. A-2, Block diagram of the modified Kalman filter. 



(see Fig. A-1) . 



Z (K) = 



Z (K) 



Z (K) 



When the processing continues we have the predicted states 
(measured state) with the covariances of the error associated 
with the predicted states (measured state) . Then by the 
likelihood ratio we can decide or H 2 * 



L 



A 







(A-4) 



Thus the decision is; 



L > 1 decide 

L < 1 decide H 2 



(A- 5) 



By the knowledge of the covariance matrix, and assumed 
Gaussian distribution of the state, L can be computed. 

The predicted measurement is determined in the 
filter procedure: 



A(K) = the covariance error of the predicted 

measurement, Y(k|k-1) (A-6) 
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Y (K) = Y(K|K-1) = 

— |J ~ 



Z (K)-’ 



(K) 

"P 



; the predicted measurement 



By defining 



r, (K) = Y (K) - t, (K) 

-1 ~p -1 



r2(K) = - t2(K) 



ti(K) ^ 






(2) 



(A-7) 



- 



nz(2) 



( 1 ) 



Thus : 



Li(K) = r^^(K) - a\k) . r^(K) i=l,2 



(Z^^^ ,Z |H^) = (2tt) ^/^exp{- I L^(K) } 



m = dimension of the measurement vector. 
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Substituting into Eq. (A-4) we get the decision law: 



L^(K) < L 2 (K) decide 

Lj_(K) > L 2 (K) decide H 2 



(A-8) 



See the block diagram in Fig. A-3 for details, 
b. The Filtering Block 

After the decision about the measured state is 
made, we have now the problem of filtering the data. The 
filtering could be straightforward if we have high confidence 
that our decision was correct, which is assumed in the 
regular Kalman filter. In our case we are not sure if our 
decision was correct, and so in the filtering process, we 
must take this into account. Because we cannot tell if our 
decision was correct or not we cannot reject or accept the 
data, we can only "weight" the data and depend on the 
computed probability that our data is correct or not. 

We have two possibilities: 

(i) A true data available, i.e., we decide 
the right decision, we will call this event D(K) - the 
data is desired, in sample K. 

(ii) A false data available, i.e., we decide 
the wrong decision, we will call this event F(K) - the data 
is ^alse, in sample K. 

So we can write: 
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f(X(K)iZ(K)) =a(K)-f[X(K) D (K) , Z (K) ] + [1-a (K) ] f (X (K) |F (K) , Z (K) ) 

(A-9) 



where 



a(K) = Prob{decision K was desired} 
= Prob{D(K) 1 Z (K) } 



(A-10) 



Assume that L 2 (K) . Then: 

a(K) = P{D(K) |Z(K) } = P{Hj_(K) |Z(K) } 
f {Z(K) 1 h^(K) }*P(Hjl(K)) 

fiZ(K) IH^(K) }-P(H^(K)) +f (Z(K) !H2(K) -P1H2(K) } 

(A-11) 



Observe that: 



P{H^(K)} = P{H2(K)} = I 



Then: 



a (K) 



1 

1 + B^(K) 



i = 1,2 



(A-12) 



where, in case of L 2 (K) > L^(K): 

^ f (Z(K) IH 2 (K) } 
^ f {Z (K) IH^ (K) } 



-( L 2 - L ^)/2 

e 



and in case of L^^(K) 



> L 2 (K), then; 



^ f (Z (K) (K) } 

^2^^^ ^ f {Z (K) IH 2 (K) ) 



-(L^-L2)/2 

e 
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Then we can conclude that; 



a(K) 



1 

1 + B(K) 



(A-13) 



where : 



. -|L, (K) - L.(K) 1/2 

B(K) ^ e 



(A-13a) 



From Equation (A-9) we see that: 



X(K|K) = E{X|Z{K)} = a(K) E{X(K) I D{K) ,Z (K) } 



+ [1 - a(K) ]E{X(K) |F(K) ,Z(K) } 



(A-14) 



i.e, the new estimate of the state will be weighted with 
the results of a filter when ass-uming true data, and the 
predicted estimate of the state (assuming false data) . The 
question now is, what kind of filter to use here? It seems 
reasonable to use a modified Kalman filter, because Kalman 
filter is the best linear sequential filter (when only the 
second order statistics are given) . We use here the notation 
"modified Kalman filter" because it is not a regular Kalman 
filter in the sense that the covariance matrices and the 
gain are data dependent (because the output is a function 
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of a) • Furthermore, there is additional selection block 



which does not appear in R.K.F. 

The only way not to have extremely time con- 
suming computations is to assume that the p.d.f. is 
Gaussian. The Kalman filter gives us the infomation that 
is used to evaluate an approximation of the p.d.f., assuming 
a Gaussian distribution. 



are data dependent, we must compute them on line rather 
than off line. 

Let's sxammarize the results first and then 

derive them: 

a. The estimator is characterized by: 



Because our gain and the covariance matrices 



A 

X(K|K) = X(K|K-1) + G^^(K)[c(K) - CX(K|K-1)] 



(A-15) 



X(K|K-1) = $X(K-1 I K-1) + AU(K-l) 



where ? is the input vector after the selection, 
b. The gain is given by: 




a(K) G(K) 



(A-16) 



where 



G(K) 



P(K|K-1) c'^[CP (K|K-1) + R] 



G(K) 



is the regular Kalman filter gain. 



i 



I 



c. Update covariance error matrix, and predicted 



error 



where 



covariance matrix: 

P(K+1|k) = $^(K|K) + Q (A-17) 

^(K|K) = a^(K)P(K|K) + [1 - a(K)]^P(K|K-l) 

+ 2a(K) [1 - a(K)]{[I - G^(K) C] (K-1 1 k-1) 

+ [I - Gj^(K)C]Q} 



P(K|K) = [I - G(K)C] P (K|K-1) 
P(k|k) = the regular Kalman filter. 



d. 

a(K) = I 

1 + 6(K)- 

^ - |L. (K)-L. (K) 1/2 

3(K) = e ' 2 

a^(K), L2 (K) are given by (A-7) 



(A-18) 



e. Initial conditions: 



P(Ol-l) = M 



(A-19) 



#'r 

( 



X(0 1-1) 



f. The predicted measured error covariance 
matrix is given by: 

A(K) = CP(K|K-1) c'^ + R(K) (A-20) 

(See Fig. A-4 for details.) 

Explanation and Notation 

1. The differences between the Kalman filter 
and this filter are: 

— The gains are data dependent, and so they 
must be computed online. 

The error covariance matrix is not the 
same as in the Kalman filter. 

— We do not use the measured vector Y but 
the vector Z. 

2. The Modified Gain - G,. : 

M 

In the regular Kalman filter the estimate 
equation of the state is given by: 

X^^NkIk) = X(K|K-1) + G(K) [? (K) - CX(K|K-1) ] 

(A-21) 

/V 

X(k1k-1) = <I>X(K-11k-1) + AU(K) 



where 



X 



(K) 



(k|k) is the update estimate state in R.K.F. 



r* 



L 




I 



Fig. A-3. Block diagram of the selection block. 




Fig. A-4. The block diagrame of the filtering procedor. 
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In our case, when we compute the state X(k|K) we must take 
into account both the updated states assuming that our 
measurement is correct, and the predicted states assuming 
that our measured vector is not correct (A-14,15), i.e.; 

X(K|K) = a(K)X^^^ (KjK) + [1 - a(K) ]X(K|K-1) 

(A-22) 

/N /N 

= X(k|k- 1) + a (K) G (K) (K) - CX(kIk-I)] 

Thus we see that the modified gain is = a{K)G(K). 

3. Covariance matrix of the error: 

The covariance matrix of the error is 

defined by: 



^(K|K) = E[e(K|K) •e'^(K|K) ] 



(A-23) 



And: 



e(K|K) = X(K|K) - X(K) 



/\ /V 

= a(K)X(KlK) + [1 - a(K) ]X(K|K-1) - 



Substituing e into ^ we get: 



J 
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X ? 



^(K|K) = a^(K)P(KlK) + [1 - a(K) ]^P(K|K-1) 

+ 2a(K) [i - a(K) ]E(e(K|K)e'^(KlK-l) ] 

where P(K|K), P(K|K-1) are the regular Kalman filter error 
covariance matrices. Now we have to evaluate the estimation 
of 

e (K| K) e'^(K| K-1) ; 

It can be shown that: [35,38]: 

e(K|K) = (I - (K)C]$e(K-llK-l) + (^ (K)V(K) 

- [I - (K)C] W(K-l) 

M 

e(K|K-l) = $e(K-l|K-l) - W(K-l) 

Thus : 

e (K[ K) e'^ (K| K-1) = { [I-(^ (K)C]$e (K-l| K-l)+<5j^(K)V(K)-Gj^(K)V(K) 

- [I-G^(K) ]W(K-1) }{$e(K-l|K-l) -W(K-l)} 

= [I-Gj^(K)C]$e(K-l|K-l)e^(K-ll 

+ Gj^(K) V(K)e^(K-ll K-D^"^- [I-Gj^(K)C]W(K-l)e‘^(K-ll K-D'd'^ 

- [I-G^^(K)C ]$e(K-l 1 k-1)W^(K-1)-Gj^^(K)V(K)W^(K-1) 

+ II - G (K)C]W(K-1)wT(k-1) . 

M 



Using the properties: 



E{e(K-llK-l) V^(K)} = 0 

E{e(K-l|K-l) w'^(K)} = 0 

E{W(k-l) V^(K)} = 0 

we get 

E{e(K|K)e'^(KlK-l) } = [I-G (K) C] 0^ (K-1 ] K-1) <5^ 

~ ~ M 

+ [I - Gj^(K)C]Q 

Substituting the result into the equation for '[ we get: 

•J;(k1k) = a^(K)P(KlK) + [1 - a (K) ] (K| K-1) 

+ 2a(K)[l - a (K) ] { [I-Gj^(K)C] (K-1 |k-1) <5^ 

+ [I-G (K)C]Q}. 

M 

4. The predicted measured error covariance 
matrix (A (K) ) : 

A(K) = E{[Y(k1k-1) - Y(K) ] [Y(K| K-1) - Y(K)]^} 

y(k1k-i) = CX(K|K-1) ( 
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Y(K|K-1) - Y(K) 



= C{X(K|k-1) - X(K) - V(K) 



From Eq. (A-24) we get 



A(K) = CP(K|K-1) + R(K) . 



These equations are simulated and tested as discussed 
in Chapter IV. 
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APPENDIX B 



THE MULTIVARIABLE NORMAL DISTRIBUTION 



In this dissertation the multivariable normal distribu- 
tion appears in many cases. As with the M.K.F., this 
appendix siammarizes the equations of such a process. 

Let's denote the random variables x^, ..., x^ by the 
column vector: 



X 



1 



X 



(B-1) 



X 



n J 



The covariance matrix A is defined as; 



o 



11 



o 



12 



a 



In 



a 



21 



A 



(B-2) 






E{ (X^ - Uj^) (Xj - Uj) ] } 



(B-3) 
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The matrix is symmetric because 



a . . = a . . 

ID Di 



(B-3) 



And the vector u (the vector mean of x) by: 



U = E 



^1 




I 

1— 1 
0 

1 


• 




• 


• j 




• 


• 




• 


X 




u 


L 




L n _ 



(B-5) 



It can be shown that, for the n-variate case, the joint 
density f (x) is given by: 



f(X) = (2tt) I A I ^-^^exp{- j(X-U)'^A ^ (X-U) } (B-6) 



where 



A I = det [A] 



(B-7) 



Also, we can write in matrix notation that 



A = E{ (X-U) (X-U) } 



(B-8) 



And : 



U = E{X} 
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As a matter of interest, we note that the n-variable 
normal moment generating function can be written as 



M(t) = exp{u'^t + t"^At/2} 



(B-9) 



where: 



t 



A 




(B-10) 



M(t) = E{exp(A)} 



(B-11) 



4 

And by direct differentiation we can find that: 






} = 



^a+b+. . . ^ 
3t^^ t^^ . - . 



t=0 



(B-12) 



The results are for the four moments, in four dimensional 
case : 



= E[X^X2]E[X3X^] + E [X^X3] E [X2X^ ] 



+ E[X^X^]E[X2X3] - 2E[Xj_]E[X2]E[X3E[X^] 



(B-13) 
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One of the interesting and useful properties of normal 
random variables is that they are invariant to linear 
transformations. In other words, a linear transformation 
of a normal random variable is a normal random variation 
too . 

The result is that for the transformation 

Y = A X (B-14) 



where: 




A = n X m matrix 



is 



= E[Y] = A U (B-15) 

= COVAR[Y] = AAa’’^ (B-16) 

Let's take a special case but very useful for bivariate 



normal distribution with: 



J 



0 



E[Xj_] = E[X2] = 



VAR[X^] = VAR[X2] = o‘ 



(B-17) 



then ; 



E[X^X2] = a p 



A = a‘ 



and we can write the p.d.f. by; 



fj^(X^,X 2 ) = — ^ \/l-P^ exp[- 



2t\o 



X ^+X.^+2pX,X 

— 5^ y ---- - 3 (B-18) 

2a^(l-p^) 



The p.d.f. of one variable given by 



f^(X) = — ^ - exp [- 



V 



27ro‘' 



2a 



(B-19) 



The conditional p.d.f. is: 






(B-20) 



Thus 
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1 






\j 2ira^(l-p^) 



p^X ^+X-^+2X.X 

exp [ ^ ^ LJ.] 

2a^(l-p^) 



(pX^+X2) 






2tto^ (1-p^) 



exp [ 5 5 -] (B-21) 

2a^(l-p^) 



r 

f 



j 



I 



APPENDIX C 



TRANSFORMATION OF RANDOM VARIABLES 

Consider the problem that arises when a random variable 
X is transformed to a new random variable y through a 
functional relationship y = h(x). The usual problem is, 
find the distribution of y, given the distribution of x. 

For a monotonic function, it can be shown that [22] , [1] : 

f^(Y) = fj^(X) |||| (C-1) 



Or by writing; 



X = h"^(y) 



(C-2) 



Then; 



f^{Y) = fj^[h"^(Y) ] ] 1 (C-3) 

This only applies when h ^ exists and it is continuously 
differentiable function of x. For this case, we must 
start from the basic idea which is; 

P(Y < b) = P(Set of X corresponding to Y < b) (C-4) 
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I 




For example, the set of x corresponding to Y ;< b is the 
intervals 



(-«»,x^) , (X2,X2) 



So by the definition of we can write: 



P(Y<b) = Fy(b) = /fj^(X)dX 



(C-5) 



where 



r = {all the values of x for which Y ^ b}. 

A very useful example of transformation is the square 
law transformation: 



y = 



2 



X 



(C-6) 



Then : 
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m 



b > 0 



(C-7) 



F^ih) = / fj^(X)dX ; 

-/b 



differentiating with respect to b, we have; 






f^(/b + /b) 

2/b 



Replacing b by y we get: 






f^(/Y) + fj^(- /Y) 

2/y 



(C-8) 



Now let's take the normal Gaussian density; then we get: 



fy{Y) = — [exp(-Y/2) + exp(-Y/2)] 

2/2ttY 



or 



f^(Y) = (27TY) exp[-Y/2] 



(C-9) 



Now, let's take the case of two variables: 



Yl = 

Y2 “ h2(Xj^,X2> 



(C-10) 
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It can be shown that if and related to and 

X 2 through continuously differentiable transforms then 

[ 1 ], [ 22 ]: 



f^(Y) = f^(X) I J(X|Y) 1 



(C-11) 



where : 



J (x I y ) = det 







(C-12) 



J is called the Jacobian of the transformation. Where the 
transformations are not one to one, we must start from the 
definition : 



Fy(b^,b2) = /^/f(X^,X2) dX^ dX2 (C-13) 

r = {set of x^,X 2 such that y^ £ b^ 
and y 2 £ ^2 ^ 

Now, by differentiating with respect to b^ and b 2 and 
setting “ Y]_' ^2 ~ ^2' proper p.d.f. 
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APPENDIX D 



SEPARATION OF SIGNALS BY PROBABILITY FILTER 
1 . Introduction 

The general problem of separation of two signals is 
well discussed in the literature and the most familiar 
are Kalman and Weiner filters. We, in general, call one 
signal "noise" and the second signal is called simply 
"signal", and the problem is to estimate the signal when 
we have a noisy measurement. The most discussed and 
developed method of separation is for linear models, i.e., 
in control theory we have a general discrete model; 



X(i+1) = <^X{i) + AU(i) + r W(i) 



(D-1) 



Z (i) = C X(i) + D V(i) 



When X is the state vector, Z is the measured vector, W 
and V are noises, U is known control input,(£),A/ T, C, D 
are matrices (may be a function of time) . 

A problem arises when the system is not linear. There 
are many possibilities to solve this problem, and the most 
widely used are the extended Kalman filter, maximum likeli- 
hood filter, and the Bayesian approach. Here we use a new 
method which will be called the "probability filter." The 
advantage of this filter is that the knowledge required for 
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this filter is less than for the likelihood or Kalman 
filters, i.e., the only required knowledge is at most the 
a priori p.d.f. of the measured signal (sometimes even the 
knowledge of the p.d.f. might be non accurate). We will 
show that this method is an extension of the method known 
as "method of moments." In the next paragraph we will 
describe briefly the method of moments and after that we 
will introduce the new approach. 

2 . Review of the Method of Moments 

The method of moments is simple and intuitive. It is 
practical in the sense that it leads to a computationally 
reasonable estimator. This method can be established without 
the a priori knowledge of its p.d.f., but it requires a 
conditional p.d.f. on the observation. In practice, it 
requires even less than this — only knowledge of the several 
first moments of the process. It yields an estimation that 
is not necessarily optimal in any sense. On the other hand, 
the estimate approaches the true value as the amount of 
data processed becomes infinite. So, on the one hand, it 
is sometimes slow in obtaining results because the amount 
of data required for the estimation is high, but, on the 
other hand, in many cases the result is guaranteed to estimate 
the true value. We conclude that, in general, this method 
simply estimates the true value, but it is not an efficient 
estimator. 
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Now let us introduce the basic concept of this method 
(according to Ref. [1] and [51]). The two references do 
not refer to the ambiguity problem which we have discussed 
above. Suppose that we have unknown vector parameter 



= [0 , . . . , 0 ] 

~ 1 n 



(D-2) 



and we have sampled data from Z(t); 



Zj Z(tj) f j 1^2/ ...^ m 



Z^ = [z^, ..., z^] 



(D-3) 



which depends on the vector parameter. Let's assume that 
the knowledge of the conditional moments from 1 to J of 
Z(t) is given, and exist, i.e., 

Ui “ = f^(0) , i = 1,2, ..., J (D-4) 

th 

Let's suppose now that we have J estimators of the J 
moments, so we can write: 



“l = 

(D-5) 

“i = 



1 
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Now we have J equations with n unknown parameters, 6 ^, 02 / 

6 ^. In general J ^ n because the functions are not 
linear. 

If we have the same number of equations as unknowns 
it can lead to a non-unique solution. So for uniqueness 
we can use more equations than the number of the unknowns. 
Note: The problem of ambiguity can be also addressed by 

the Modified Kalman filter (see Appendix A) . 

Example: D-1. (Refer to Fig. D-1.) 

If we are given only two equations, f^(0) and f 2 ( 6 )/ with two 
parameters, 0 ^ and we have two possibilities for the 

vector 0: 0^^^ and 0^^^. To decide what is the acceptable 

solution, we must have at least one more equation. With 
the addition of a third equation we accept one of the 
solutions. Note that the third equation does not lead to 
the same solution as the first two equations, because of 
the error in estimating the moments . Because in many of 
the cases, the estimation is an increasing function of J, 
the last equation is used only for decisions and not for 
estimation, i.e., our solution is: 



1*3 <9 


(D) 


- S 3 I 


< lf3(0^^^ 


/N 

“3 


l*3<? 


( 2 ) 

) 


- a3| 


< |f3(0^^^ 


/N 

“3 



choose 0 



( 1 ) 



(D- 6 ) 



choose 0 



( 2 ) 



The block diagram of the method of moments is given in 
Fig. D-2. 
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Fig. D-L.. The concept of determining the solution for two 
unknowns with three non-linear equations. 
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Fig. D-2. Block diagram of the method of moments estimator. 



1 



In most of the cases the computational method of 
estimating the moments are simply by the definition of 
the moments , i . e . , 



/N 




m 




i=l 



(D-7) 



It is not difficult to generate examples in which the 
method of moments does not even lead to an acceptable 
estimator (see Ref. [1]). In the following section we will 
see a new method which is a generalization of the method 
of moments . 

3 . The "Probability Filter" 

In the method of moments we saw that when one wants 

t h 

to estimate the parameters the first moments of the 

observed data is taken and we equate these to the J functions, 
the moments, and solve for the unknown parameters. Let's 
now suppose that J -> <». We choose the parameter in such 
a way that these infinite equations will lead to minimum 
error in some sense. For example; 

00 

\ [a^ - f^(9)]^ • Wj^ min (D-8) 

i=l 

when ^ 0 . In general, w^ is chosen as a decreasing 
function of i, for example w^ = l/var{a^ - f^(9)}. The 
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method of moments is one example of D-8 , i.e., by choosing: 



w . 

1 



1 



i = 1 



. . . , n 



w . 

1 




1 



n / « • • / 



J 



(D-9) 



0 



i > J 



We get exactly the method of moments discussed before. Now, 
let's understand what we mean by Equation D-8. If we take 
> 0 V- i, then we are trying in some sense to minimize 
the difference between the conditional p.d.f. and the 
measured p.d.f. That is, we can modify Equation D-8 by: 



when ^ 2 ^^^ denotes the estimation of the p.d.f. of the data. 
(We will see later that for the probability filter we do 
not even have to estimate the p.d.f.) In the last equation 



we try to choose 9 in such a way that the norm of the differ- 
ence between the two p.d.f. 's is minimized. 



I |f2(Z|0) - fj(Z) I I 



^ min 



(D-10) 



where 




f2(z(e) ) 
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Example D-2. (Refer to Fig. D-3.) 

Let's asstime that we have only one parameter to esti- 
mate; from the figure we choose 0^^ and not 0 2 because 0^^ 
gives a theoretical p.d.f. which is closer to the measured 
p . d . f . than 0 2 • 

Let us now choose the norm: 

00 

1=1 / [f2(Z|0) - dZ ^ min (D-11) 

— 00 

That is, we choose 0 in such a way that it minimizes the 
difference between the two p.d.f. 's in a mean square sense. 

Now, let us take the derivative of I with respect to 
0, and set it equal to zero to find the minimum, i.e., a 
necessary condition for minimization is; 



0 = 



31 

/N 

30 . 
1 



- - 3f (Z|0) 

/ [f„(Z|0) - f„(Z)] * ^ dZ 

^ ^ 30. 

5 1 



31 

/\ 

30. 

1 



3f„(Z|0) . “ 3f„(Z|0) .. 

/ ^ fyiZ\Q) dZ - / ^ f,(Z) dZ 



30 



30, 



gi(0) - 



Bf^CZI©) 



30 . 
1 



= g.(0) - 

J. ^ 



3f2(Z|0) 

/A 

30 . 

1 



(D-12a) 
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1 



where the symbol means time average. The time average 

can be approximated by; 

, T 3f7(Z(t) |6) 

^ dt (D-12b) 

^ 0 30j_ 

In our problem we try to estimate the amplitudes 
received at each receiver (which is porportional to the 
RCS of the target, see Chapter IV. A) . Since we restrict 
ourselves to surface targets (which means low dynamic) , 
the RCS of each target (ship and chaff) can be assumed to 
be stationary compared to the missile's flight time. 

Hence we can use "T" , the time interval at the last 
equation (D-12b) , small enough compared to the missile's 
flight time. 

Eq. (D-12,D-12a) follow because; 

3f„(Zl0) 

^ fy{Z\e) dZ 

30 . z ~ 

1 

is a known function of 0 and can be computed a priori . 

Also 



9i<?) = / 

— oo 




/ f (Z) dZ 

3«i 



3f,(Z|6) 

is the mean value of the function with respect to 

30 . 
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Z. Because we are assuming ergodic processes, the ensemble 
mean and the time mean are the same . So we can write : 

[Z(t)|e] 

g (0) = 0 (D-13) 

^ 39 . 

1 

Because we have n parameters, we have n equations. 

Note that also in this case as in the method of moments , 
the equations are not linear and we may get a non-unique solu- 
tion. Thus we may have to add more equations to the set of 
equations D-13; for example, the equations of moments as 
discussed previously. Alternately we propose to use the 
M.K.F. (Appendix A) with all the possible solutions. 

A block diagram of the probability filter is given in 
Fig. D-4. The input signal Z(t) is substituted into 

3f2(Z(t) 19 ) 



and the function is averaged by the block following. Also, 

<A. 

g(9) is produced. If the difference between the two is not 

/N 

zero, the value of 9 is changed until the difference between 
the two approaches zero. The initial condition for the 
integrator can be determined by the method of moments or 
as in the Kalman filter, the statistical mean of the 
parameters . 
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ej> 




Fig. D-4. Block diagram of the probability filter. 



Example D-3 » 

In this example we introduce a simple case of estimating 
a parameter by the probability filter. We choose the 
exponential distribution. The p.d.f. is given by: 



f 2 (Z 1 0) = 6 e"®^ Z ^ 0 


(D-14) 


We can get the equation for the probability filter; 

-ez „ . -ez 

36 ■ ® ® 


(D-15) 


" 3f,(Z|0) 

9i(e) = / — Se f2(z|9) dz = + T 

0 

Thus gj^(0) has a very simple expression, i.e.. 


(D-16) 


g^(9) = + 0.25 

Thus the equation for the probability filter is, from 
Eq. (D-13) , 


(D-17) 


N 

, -0Z. ^ 

+0.25 + ^ 1 e ^(9Z^ - 1) = 0 

i=l 


(D-18) 


The last equation is the probability filter equation. 


It 



is a non-linear equation and we must solve it numerically 

/N 

for 9 . 
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Until now we have seen the direct implementation of 



the filter, i.e., determine the probability filter equation 
and solve for the parameters. Direct implementation of 
the filter is only one possible way of implementation. 

There are other possible implementations which give the 
filter its uniqueness. We introduce two possibilities of 
using the filter. In many cases we do not have the exact 
p.d.f. of the process but we have an approximation of it 
and a reference signal which represents the signal (e.g. , 
by recording the signal previously) . Wg give here two 
ways of using the probability filter when the exact p.d.f. 
is not known. 

(i) The solution when only the approximate p.d.f. is 
given and a reference signal is available. 

This filter gives us the opportunity to estimate parameters 
even if our information about the data is not precisely 
given; that is, if we have a record of this kind of data, 
we can compare the two systems very easily. 

Let's return to Eq. D-12b. We have: 



oo 



(Z 1 0) . 




/ 



af^czje) 



•f2(zl 6)dZ 



/ — <3Z 

80 . ^ 



— OO 




set 



0 



(D-19) 
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Let us assume that f^{Z\Q) is only an approximation of the 
p.d.f., so we cannot "trust" this given p.d.f. completely. 
Like in the Kalman filter, we assume an approximate model 
of the system, and by closing the loop we guarantee that 
the estimator behaves quite well even if our model does 
not precisely describe the real model. 

Equation D-19 gives us the possibility to close the 
loop as in a Kalman filter. Eq. D-19 has two terms. The 

/N 

first one is computed "off" line if the p.d.f. f 2 (Zj 0 ) is 
completely known. If the p.d.f. is unknown exactly, then 
the technique of comparing two terms is the best that we 
can do. We see that the two terms are almost the same. 

Assume we have a reference signal X with the p.d.f.: 



/s 



f^(Xl0) 



= f2(Xl0) 



(D-19a) 



Thus : 



9f„(Zl0) ^ 9f„(Xl0) 

/ f„(z|6)dz = — 

90 . ^ ~ 90 _ 

’ 1 1 



Thus we can write: 





9f2(Xl0) 





(D-19b) 
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X 



= x(t,0) 



Z = Z(t,0) 



The first terro comes from the reference system, so by 
computing the two terms we can assume that when the 
difference between the two terms becomes zero, 0 is close 
to the true 0 . 

Figure D-5 demonstrates the idea. We measure the data, 

/N 

Z(t,0), transform the data to get the function h^ ( Z , 0 ), where 



/V 

h^( w , 0 ) 




(D-20) 



and g(oj|0) is an approximation of f(a)|0) that has been 
chosen for its analytic form. 

In the original form of the probability filter equation 
(D-12a) , we get an error when f (y| 0) is not known exactly. 
Thus by implementing the probability filter equation (D-12a) 
we get a biased estimator. On the other hand if we have 
a reference signal, in case that is not known exactly, 
we get: 







h^(X(r,0),0) - h^ (Z (t,0) ,0) = e^ 



(D-21a) 
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Figure D-5. Probability filter when a reference is given. 




1 



For this case when 6=0: 




(D-21b) 



oo 



oo 




— 00 



— CO 



0 



This result follows because of the assumption that we have 
a reference signal with the same p.d.f. as the signal, 

(D-19a) . 

Thus where the reference signal has the same functional 
equation as the measured signal, we get an unbiased estimator, 
independent of hj^(C/9)/ which is a function of the assumed 
p.d.f. However the feedback process can only be expected 
to converge when g(wl6) is a reasonable approximation to 



In the regular "probability filter" we get a bias 
estimator when the p.d.f. is not known exactly. 

(ii) Indication of Change in a System 
Sometimes we like to check if parameters in a system have 
been changed. By taking a reference signal (from similar 
systems with nominal parameters or by recording the output 
from the system from the past) , we can easily check if some 
parameters have been changed. 



f2(<^l6) • 
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i 



Indication of change in the system is that: 



le^l > (D-23) 

where is a constant dependent upon the acceptable false 
alarm. The block diagram is given in Fig. D-6 . 

By adding more calculations we might determine the 
value of the parameters that were changed, by simply 
applying the probability filter. 

The probability filter implementation for the unresolved 
targets and the simulation of the probability filter are 
given in Chapter V. 
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APPENDIX E 



THE MAXIMUM LIKELIHOOD ESTIMATOR 



A reasonable estimate of a parameter is one which 
makes a given observation most likely, i.e., choose for 
the parameters 0 , that value which maximizes the conditional 
p.d.f. f(z|0). The p.d.f. f(z|0) is called the likelihood 
function. Because of the exponential character of many 
density functions, it is convenient to deal with the natural 
logarithm of the likelihood function - ln{f(z|0)}. Clearly, 
because f(z|0) ^ 0, In f{(z|0)} exists, so by choosing 0 to 
maximize ln{f(z|0)} it maximized also f(z|0). A necessary 
condition that 0 maximize the likelihood function is: 




In f (Z| 0) 



0 i=l,2,...,m (E-1) 



By setting all the m equations with m unknowns, one can 
solve for the unknown vector 0. It can be shown that in 
many cases the maximum likelihood function is a biased 
estimator, but the bias tends to zero [30]. Also the 
variance of the error tends to zero as numbers of samples 
increase. We can say that the likelihood estimator is 
asymptotically efficient. 
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APPENDIX F 



ESTIMATORS FOR SYSTEM DESCRIBED BY 
DIFFERENTIAL/DIFFERENCE EQUATION 

1 . INTRODUCTION 

In this dissertation we deal with estimators. The 
purpose of this appendix is to summarize the approaches 
given in the literature and establish the symbols and 
terminology used. There are two major approaches: the 

discrete and the continuous time. For our case, we deal 
with discrete time process. 

The difference equations for the discrete time and the 
measurement is given by: 



X(K+1) = (X (K) ,K+1,K) + r (X(K) ,K) W(K) ) 



(F-1) 



Y(K) = h(X(K) ,K) + V(K) 



where 



ijj,r,h - are known matrices 



W(K), V(K) - are independent sequences of indepen- 
dent random variables (not Gaussian 
in general) 



The only known case when the equations can be effectively 
solved (computed) exactly is the linear process plus 
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Gaussian noise. This solution was obtained by Kalman [2] 
and Kalman and Buoy [3] . The known existence of this exact 
solution led to the approximate solution for the case where 
the plant and the measurements equation are almost linear 
and the noise processes are nearly Gaussian [10] . But the 
Gaussian approximation and/or linearization are not adequate 
for many problems due to strong nonlinearity, as well as for 
non-Gaussian noise process. This has led to the development 
of many other kinds of approximations, for example, Gaussian 
sum approximations [5] , [6] . A good reference for many 

approximations can be found in the proceedings of the 6th 
nonlinear estimation theory symposium [7]. 

From equation (F-1) we see that the frequent problem 

is : 

The measurement noise and the plant noise are 
additive, with at most a state dependent multiplier, when 
the noise is in the form of Brownian motion. 

2. LINEAR SYSTEM - KALMAN ESTIMATOR 

Let's assume that our system described by a linear 
difference equation which is given by [31]: 



X(K+1) = $X(K) + AU(K) + W(K) 



(F-2) 



Y(K) = CX(K) + V(K) 



where 
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X(K) 

U(K) 

W(K) 

V(K) 

Y(K) 



are the states of the system at time K 

the control input to the system at time 
K (it is a known, deterministic) 

the random forcing input 

the measurement noise 

the measured vector 



The statistical properties of the variables are given by: 

(1) The measurement noise: 



E[V(K) ] = 0 K ^ 0 

E[V(K)V^(j) ] = R6^j k, j ^ 0 

(2) The initial state: 

E[X(0)] = X^ 

COVAR [X(0)1 = M 

(3) The forcing random process: 

E[W(K)] = 0 

E[W(K)w'^( j) ] = k, j > 0 



(F-3) 



(F-4) 



(F-5) 
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I 



( 4 ) The 


correlation between the variables: 


(a) 


The measurement noise and the initial states 
are uncorrelated: 

E[X(0)v'^(K) ] =0 K>0 (F-6) 


(b) 


The random forcing input and the initial 
state are uncorrelated: 

E[X(0)w'^(K) ] =0 K^O (F-7) 


(c) 


The forcing random process and the measurement 
noise are uncorrelated: 

E[W(K)v'^(j) ] = 0 K,j >0 (F-8) 



The requirements for the filter are: 



(a) 


The estimator is characterized by linear 
relationship. 


(b) 


Unbiased estimator. 


(c) 


The estimator minimizes the trace of the 
error covariance matrix, i.e. : 

Trace {P(K K) } min (F-9) 


where 
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« 

I 

I 



iJ 



i 



P(K|K) = C0VAR{X(K|K) -X(K)} 

/N 

X(k[k) is the estimated state depending on the 
input data up to and including sample K. 



The results of the Kalman estimator are given by [ 31 ] : 



X(K|K) = X(K|K-1) + G(K) [Y(K) - CX(K|K-1)] 

/N 

X(K|K-1) = $X(K-11k-1) +Au(K-1) (F-10) 

X(Ol-l) = X(0) 



where : 



G(K) = 


P (K 1 K-1) [CP (K 1 K-1) c'^+R] 


p(k1k) 


= [ I - G(K)C] 


P(K1k-1) 


P (K+1 |K) 


= $p(k1 k)a’^ 


+ Q 


P(Ol-l) 


= M 





Notes ; 

(i) The Kalman filter is the optimum linear estimator 
for the given knowledge of the second order 
statistics only. 



217 



(ii) There are other estimators, like nonlinear, which 
are superior to the Kalman filter for the general 
case . 

(iii) If all the variables are Gaussian, the optimal 
minimum variance estimator turns out to be the 
Kalman filter. 

3. THE BAYESIAN APPROACH 

The Bayesian approach simply takes the expected value 
of the state, given all the measured data, i.e. ; 



X(k|K) = E{X(K) I Y(0) , 



• • • / 



Y(K) } 



(F-11) 



If we define: 



Y 



(K) A 



[Y(0) , 



Y(K) ] 



• • • / 



Then (F-11) can be formed as: 



X(K|K) = E[X(K) [Y^^^ 



( K ) 

If we have the p.d.f. of the state given Y 



, we can 



find the expected value of X, i.e. : 



X(k1k) = / X(K) f[X(K)|Y^^h dX(K) 

[X] 



(F-12) 
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So the problem is to find the p.d.f. 



f [X(K) I ] = 



f [X(K) ] 



f [Y(K) ] 



f[X(K) I f [Y(K) |X(K)] 

f [Y(K) I ] 



= / f [X(K-D 1 f [X(K) |X(K-1) ]dX(K-l) 

[X] 

= / f [X(K) f [Y(K) lX(K) ]dX(K) 



[X] 

f [X(0) ] = f fX (0) ] 



Those equations, theoretically, give the complete solution 
of the estimation. The problem is that the integration 
can not be carried out in closed form. Thus, in most of 
the practical cases this solution is not used. 
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